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INTRODUCTION  AND  ACKNOWLEDGEMENTS 


Welcome!  This  is  the  second  Alberta  workshop  on  applied  statistical  and  biometrical  applications.  We 
are  again  especially  impressed  by  the  wide  response  to  our  invitation  and  the  diversity  of  interests.  Your 
active  participation  will  make  this  year's  workshop  exceed  all  expectations,  when  even  compared  to  last 
year's  success. 

We  are  fortunate  to  have  received  the  continued  support  of  the  Alberta  Agricultural  Research  Institute 
(AARI)  and  the  AARI  Board.  Dr.  Ralph  Christian,  Ms.  Emily  Samoil,  and  all  AARI  staff  have  been 
willing  to  assist  us  in  every  way. 

Special  debt  is  owed  many  persons  who  contributed  time  and  effort,  often  having  to  work  within  tight 
schedules.  Mr.  Phil  Henry,  Alberta  Environmental  Centre  (AEC),  once  again  has  been  responsible  for  the 
computer  graphics  and  layouts,  registrations  and  generally  keeping  the  process  in  control.  Mr.  John  Kirtz, 
Dr.  Alan  Tong,  Mr.  M.  Liu  and  Mr.  Serge  Dupuis  kindly  agreed  to  assist  with  procuring  and 
demonstrating  hardware  and  software.  Dr.  Alan  Tong,  and  his  staff,  are  extended  special  thanks  for 
assisting  in  the  compilation  of  papers  for  the  proceedings;  Ms.  Grovenore  Flato,  AEC,  contributed  valuable 
help  with  wordprocessing  and  preparation  of  manuscripts.  Our  sincere  appreciation  to  Ms.  Janet  Smalley 
and  the  support  staff  in  the  Beef  Cattle  and  Sheep  Branch  of  Alberta  Agriculture,  those  in  Animal 
Sciences  Division,  AEC,  and  Ms.  Gloria  Horon,  Central  Duplicating,  AEC,  who  so  kindly  met  our 
deadlines  on  several  occasions. 

It  has  been  very  gratifying  to  have  received  the  valuable  support  and  participation  from  those  contributing 
presentations.  Through  their  willingness  to  share  informative  topics  of  interest,  and  program  planning 
input  from  Milton  Weiss,  Alan  Tong  and  Len  Kryzanowski,  we  believe  an  especially  valuable  workshop 
program  has  again  come  together.  Join  us  in  extending  an  "Alberta  Welcome"  to  our  visiting  participants, 
Drs.  George  Milliken,  Kansas  State  University  and  Richard  Gilbert,  Battele,  Northwest  Laboratories, 
Richland,  Washington. 

We  are  particularly  aware  of,  and  acknowledge,  the  support  we  have  received  from  our  respective 
Departments,  the  Executive  Director  (AEC),  Division  Directors  and  Branch  heads.  To  those  we  may  have 
mistakenly  overlooked,  please  accept  our  apologies. 

Whether  there  should  be  another  workshop  like  this,  and  last  year's,  will  be  a  topic  for  discussion  on 
Tuesday  morning  (October  22,  0900  h).  Regardless  what  consensus  may  result,  we  have  been  privileged 
to  be  associated  with  many  fine  people  and  make  new  friends,  whom  we  might  never  have  had  the 
opportunity  of  knowing. 

L.  Zack  Florence 
Animal  Sciences  Division 
Biometrics  Section 
Alberta  Environmental  Centre 

Laki  Goonewardene 
Animal  Industry  Division 
Beef  Cattle  and  Sheep  Branch 
Alberta  Agriculture 
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The  Relationship  Between  Methods  of  Competing  Risks  Analysis 
and  Methods  of  Survival  Analysis,  With  an  Application  to 
Entomology 

G.  B.  Schaalje 
Research  Station,  Agriculture  Canada 
Lethbridge,  Alberta  TIJ  4B1 

Introduction 

Time-to-failure  data  arise  in  many  fields  of  application  including  medicine  (survival  times  and 
remission  times  of  patients  receiving  various  treatments),  industry  (lifetime  testing  of  products  and 
components),  plant  science  (days  to  germination  or  heading  of  various  varieties),  and  entomology 
(development  times  under  different  temperature  regimes).  'Time  to  failure'  is  used  here  as  a  generic  term 
because  a  'failure'  could  be  the  achievement  of  any  kind  of  endpoint,  and  'time'  could  be  any  kind  of 
non-negative  quantity  (eg.  distance).  A  unique  aspect  of  such  data  is  that  some  of  the  observations  may 
be  'censored'.  That  is,  for  some  reason  not  related  to  the  treatment,  an  individual  may  not  reach  the 
failure  endpoint  before  either  the  study  is  terminated  or  the  individual  is  lost  to  followup.  The  time- 
to-censoring  data  still  contain  information  about  the  failure  process,  and  a  characteristic  of  appropriate 
analysis  methods  is  that  they  account  for  this  information  in  addition  to  the  information  contained  in  the 
noncensored  data. 

This  paper  primarily  discusses  the  analysis  of  time-to-failure  data  when  there  are  two  or  more 
competing  causes  of  failure,  usually  referred  to  as  competing  risks  analysis.  Methods  for  analysing 
time-to-failure  data  when  there  is  only  one  cause  of  failure  (survival  analysis)  are  reviewed,  and  it  is 
demonstrated  that  these  methods  are  surprisingly  difficult  to  extend  to  the  case  of  multiple  causes  of 
failure.  The  cautions  that  must  be  observed  when  analysing  data  with  two  or  more  causes  of  failure  are 
discussed,  and  an  example  from  entomology  is  presented.  Software  for  both  survival  analysis  and 
competing  risks  analysis  is  briefly  discussed  as  well. 

Survival  Analysis  for  One  Failure  Type 

The  goals  of  survival  analysis  are  to  characterize  the  failure  process  in  useful  ways,  and  to  infer 
the  effects  of  treatments  and  other  covariates  on  failure.  The  various  methods  of  survival  analysis  are 
centered  around  the  survivor  function  [S(t)],  the  density  fimction  [f(t)],  or  the  hazard  function  [h(t)] 
associated  with  the  failure  times. 
Survivor  Function  Methods 

The  S(t)  gives  the  probability  that,  for  a  specified  time,  an  individual's  time  to  failure  (survival 
time)  is  at  least  that  long.  Methods  for  estimating  S(t)  for  a  homogeneous  group  of  individuals  are 
intuitive  and  simple,  the  nonparametric  product-limit  (Kaplan-Meier)  and  actuarial  estimators  being  the 
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most  common.  Statistics  such  as  the  logrank  statistic  are  available  for  testing  the  equality  of  S(t)  for  two 
or  more  groups.  As  long  as  there  are  only  a  few  discrete-valued  covariates,  stratification  of  the  data  can 
be  used  to  adjust  for  covariates.  The  SAS  program  LEFETEST  and  the  BMDP  program  IL  are  useful  in 
estimating  S(t)  and  computing  related  statistics. 
Densitv  Function  Methods 

The  use  of  f(t)  as  the  basis  for  data  analysis  is  common  when  there  is  reason  to  believe  that  a 
parametric  distribution  such  as  the  WeibuU  or  Lognormal  is  a  good  model  for  the  failure  times.  The 
advantage  is  that  the  data  can  be  summarized  in  a  few  fitted  parameter  values.  Maximum  likelihood 
methods  are  used  to  fit  these  models  to  failure  time  data,  and  two  or  more  treatment  groups  can  be 
compared  using  likelihood  ratio  statistics.  Another  advantage  is  that  regression  models  for  failure  time 
data  with  discrete  or  continuous  covariates,  usually  called  'accelerated  failure  time  models',  can  be 
formulated  in  terms  of  parametric  density  functions.  Although  the  documentation  is  a  bit  confusing,  the 
SAS  program  LIFEREG  and  the  BMDP  program  2L  can  be  used  to  carry  out  all  of  these  analyses. 
Hazard  Function  Methods 

The  h(t)  gives  the  instantaneous  rate  of  failure  at  a  specified  time  given  that  the  individual  has  not 
failed  prior  to  that  time,  and  is  related  to  the  risk  of  failure  to  which  an  individual  is  subject. 
Nonparametric  methods  can  be  used  to  estimated  h(t)  for  a  homogeneous  set  of  individuals.  A  regression 
method  (proportional  hazards  or  'Cox'  regression)  based  on  h(t)  which  combines  some  of  the  advantages 
of  both  nonparametric  and  parametric  analysis  has  been  developed.  Under  the  proportional  hazards  model, 
it  is  assumed  that  the  hazard  functions  for  all  individuals  are  proportional,  with  the  proportionality  constant 
determined  by  a  vector  of  covariates  x.  The  most  common  model  for  the  covariates  is  the  loglinear  model 
given  by 

h(t;x)  =  ho(t)  exp(b'x) 

where  ho(t)  is  the  'baseline'  hazard  function  for  individuals  with  x  =  0.  The  convenient  thing  about  this 
model  is  that  ho(t)  contains  no  information  about  b,  and  is  therefore  arbitrary.  Thus  conditional  maximum 
likelihood  can  be  used  to  estimate  b  without  any  restrictions  being  placed  on  the  form  of  h^Ct).  A  unique 
and  useful  aspect  of  this  model  is  that  time-varying  covariates  can  be  easily  incoiporated.  Also, 
stratification  can  be  used  to  accomodate  nonproportional  hazards  within  the  Cox  regression  framework  in 
some  cases.  The  SAS  program  Llhh  l  hST  and  the  BMDP  program  IL  can  be  used  to  obtain  a  nonsmooth 
estimate  of  h(t),  and  the  IMSL  subroutine  HAZST  can  be  used  to  obtain  a  smooth  estimate  of  h(t). 
Proportional  hazards  regression  can  be  carried  out  with  the  BMDP  program  2L  and  the  (soon  to  be 
available)  SAS  program  PHREG. 
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Relationships  Among  the  Functions 

All  of  these  analyses  are  equivalent  in  the  sense  any  of  the  three  functions  can  be  detemiined  from 
any  other  of  the  three  functions  upon  which  the  various  survival  analyses  are  based.  Some  of  the 
equations  relating  the  functions  are: 

S(t)  =  1  J  f(s)ds 
S(t)  =  exp[  J  h(s)ds] 
h(t)  =  f(t)/[l  J  f(s)ds] 
f(t)  =  h(t)  exp[-J'  h(s)ds] 
f(t)  =  -dS(t)/dt 

The  choice  of  an  analysis  method  depends  on  how  much  is  known  about  the  distribution  of  failure 
times,  on  the  nature  of  the  covariates,  and  on  the  goal  of  the  analysis. 

Competing  Risks  Analysis  for  Multiple  Failure  Types 

The  goals  of  competing  risks  analysis  are  to  infer  the  effects  of  treatments  and  covariates  on 
specific  causes  of  failure,  to  infer  relationships  among  the  failure  types,  and  to  estimate  failure  rates  due 
to  the  remaining  causes  if  some  of  the  causes  were  to  be  removed.  Methods  of  competing  risks  analysis 
to  approach  these  goals  have  been  developed  based  on  various  generalizations  of  the  density  function  and 
the  hazard  function.  Generalizations  of  the  survivor  function  also  play  a  role  in  some  competing  risks 
analyses,  but  usually  as  a  secondary  quantity.  Hence,  these  generalizations  wiU  not  be  discussed  here. 
Two  approaches  based  on  different  generalizations  of  the  density  function,  the  latent  failure  time  model 
and  the  mixture  model,  will  be  discussed  along  with  an  approach  based  on  generalizations  of  the  hazard 
function  known  as  cause-specific  hazard  functions.  In  this  discussion,  attention  will  be  restricted  to  the 
case  of  two  competing  causes  of  failure. 
Latent  Failure  Time  Model 

This  is  the  oldest  approach  to  analysis  of  competing  risks,  and  in  some  ways  would  seem  to  be 
the  most  natural  generalization  of  the  density  function  approach  for  one  cause  of  failure.  Under  this  model 
it  is  assumed  that  each  individual  possesses  latent  (or  potential)  failure  times  (R  and  S)  for  the  two  causes, 
with  joint  density  function  f(r,s).  For  convenience  it  is  also  assumed  that  it  is  impossible  to 
simultaneously  fail  due  to  both  causes.  Since  the  occurance  of  a  failure  due  to  either  of  the  causes 
precludes  the  other,  the  actual  cause  of  failure  is  determined  by  that  with  the  smallest  latent  failure  time. 
The  observable  failure  time  T  is  then  given  by 

T  =  min(R,S). 

If  it  can  be  safely  assumed  that  the  two  latent  failure  times  are  independent  with  density  functions  fR(r) 
and  fs(s),  then 

f(r,s)  =  f^Cs)  fs(s). 
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Various  parametric  forms  have  been  proposed  for  f(r,s),  and  estimation  methods,  regression  models,  and 
methods  for  evaluating  what  would  happen  if  some  of  the  causes  were  to  be  removed  have  been 
developed.  There  is  a  major  unresolvable  problem  with  this  approach,  however.  Because  the  latent  failure 
times  are  not  aU  observable,  f(r,s)  cannot  be  estimated  without  making  unverifable  assumptions.  Any 
inferences  that  are  made  using  the  latent  failure  time  model  also  suffer  from  this  difficulty.  This  whole 
problem  is  sometimes  referred  to  as  the  'non-identifiability  problem'  of  competing  risks  analysis. 
As  a  simple  example,  consider  the  following  hypothetical  data: 


no.  dying  no.  dying 


time 

no.  living 

(cause  R) 

(cause  S) 

0 

1000 

0 

0 

0.5 

600 

400 

1 

540 

60 

1.5 

270 

270 

2 

162 

108 

2.5 

6696 

3 

13 

53 

3.5 

0 

13 

Assuming  all  events  to  occur  at  discrete  points  in  time,  a  discrete  joint  density  function  which  fits  these 
data  perfectly  is: 


S 


R 

1 

2 

3 

4 

tot 

0.5 

.040 

.144 

.173 

.043 

.400 

1.5 

.030 

.108 

.130 

.032 

.300 

2.5 

.018 

.064 

.077 

.019 

.178 

3.5 

.012 

.044 

.053 

.013 

.122 

tot 

.100 

.360 

.433 

.107 

1.0 

However,  another  joint  density  function  which  fits  these  data  perfectly  is: 
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s 


R 

1 

2 

3 

4 

tot 

0.5 

.010 

.010 

.020 

.360 

.400 

1.5 

.010 

.010 

.010 

.250 

.280 

2.5 

.010 

.008 

.004 

.092 

.114 

3.5 

.040 

.100 

.053 

.013 

.206 

tot 

.070 

.128 

.087 

.715 

1.0 

There  are,  in  fact,  an  infinite  number  of  joint  density  functions  which  fit  these  data  perfectly.  If  some 
additional  (but  unverifiable)  assumptions  about  the  form  of  the  distribution  of  failure  times  can  be  made, 
a  unique  best-fitting  distribution  can  sometimes  be  found.  For  example,  if  the  independence  assumption 
is  reasonable,  the  first  joint  density  function  is  the  unique  best-fitting  density  to  these  hypothetical  data. 

My  personal  feeling  is  that  the  latent  failure  time  model  should  only  be  used  when  there  is  good 
reason  to  believe  the  independence  assumption,  such  as  in  industrial  reliability  testing  of  equipment 
composed  of  several  independent  components.  Such  independence  seems  unliklely  in  most  biological 
applications. 

If  independence  is  assumed,  the  SAS  program  LIFEREG  and  the  BMDP  program  2L  can  be  used 
to  fit  the  distributions  and  carry  out  regression  analyses.  To  use  these  programs  to  analyse  the  latent 
failure  times  for  one  cause,  failure  times  due  to  the  other  cause  must  be  treated  as  censored  data. 
Mixture  Model 

In  this  model,  the  failure  types  and  times  are  expressed  in  terms  of  a  stochastic  mechanism  which 
determines  the  failure  type  from  the  outset,  and  conditional  distributions  of  failure  times  given  the  type 
of  failure.  The  name  comes  from  the  fact  that  the  unconditional  distribution  of  failure  times  is  a  mixture 
of  the  conditional  distributions.  For  two  failure  types,  a  binomial  distribution  [B(n,p)]  serves  as  a  model 
for  the  mechanism  of  choosing  a  failure  type,  and  density  functions  [gR(r)  and  gs(s)]  serve  as  models  for 
the  conditional  distributions  of  failure  times.  This  model  does  not  suffer  from  nonidentifiability  problems, 
and  regression  models  for  evaluating  the  effects  of  treatments  on  the  failure  type  mechanism  and  the 
conditional  failure  time  distributions  have  been  developed. 

The  greatest  weakness  of  the  mixture  model  is  that  although  it  is  identifiable,  it  may  be 
misleading.  Consider,  for  example,  an  industrial  application  in  which  a  piece  of  equipment  with  two 
independent  critical  components  is  being  tested,  so  that  the  independent  latent  failure  time  model  happens 
to  be  valid.  Suppose  that  a  treatment  is  applied  which  affects  only  one  of  the  components.  It  would  be 
desirable  if  a  statistical  analysis  would  indicate  that  only  the  one  component  was  being  affected  by  the 
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treatment;  the  latent  failure  time  model  as  applied  to  these  data  would  show  just  that.  The  mixture  model, 
however,  would  very  likely  indicate  that  conditional  distributions  of  failure  times  due  to  both  of  the 
components  are  affected  by  the  treatment  because  in  fact  both  conditional  distributions  would  be  affected 
in  such  a  case.  The  table  below  gives  the  true  means  of  the  distributions  of  latent  failure  times  and  of 
the  conditional  distributions  of  observed  failure  times  for  three  specified  independent  latent  failure  time 
processes.  Only  the  distribution  of  latent  failure  times  for  failure  type  2  varied  among  examples,  but  the 
means  of  the  conditional  distributions  for  both  failure  types  varied  among  examples. 


Example 

Failure  Type 

Mean  of 
Latent  Failure  Times 

Mean  of 
Conditional  Failure  Times 

1 

1 

27.2 

27.1 

2 

36.8 

30.5 

2 

1 

27.2 

25.6 

2 

27.2 

25.6 

3 

1 

27.2 

24.7 

2 

24.7 

24.0 

If  there  were  no  censoring,  the  model  could  be  fitted  in  a  straightforward  way.  The  binomial  p  would  be 
estimated  by  the  proportion  of  individuals  subject  to  failure  type  1 ,  and  the  observed  failure  times  for  each 
failure  type  would  be  used  in  ordinary  ways  to  estimate  gR(r)  and  gs(s).  If  there  is  censoring,  the  EM 
algorithm  must  be  used  for  estimation  and  a  custom  computer  program  is  required. 
Cause-Specific  Hazard  Model 

This  model  is  based  on  instantaneous  failure  rate  functions  [hi(t)  and  h2(t)]  for  each  of  the  failure 
types.  These  functions  can  be  estimated  and  analysed  in  exactly  the  same  way  as  ordinary  hazard 
functions  using  the  same  software.  The  only  modification  necessary  is  that  failure  times  for  failures  of 
types  other  than  the  one  under  consideration  must  be  regarded  as  censored  observations.  This  analysis 
is  not  subject  to  nonidentifiability  nor  misinterpretation  problems,  so  long  as  it  is  kept  in  mind  that 
conclusions  based  on  these  analyses  apply  only  under  the  prevailing  study  conditions,  with  all  failure  types 
operative.  A  further  benefit  is  that  relationships  among  the  failure  types  can  be  studied  via  the  use  of 
time -varying  covariates  if  'risk-indicator'  variables  for  some  of  the  failure  types  are  available. 
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Relationships  Among  the  Models 

Relationships  among  the  various  competing  risks  models  are  discussed  by  Gail  (1975),  Gail 
(1982),  and  Prentice  et  al.  (1978). 

Example 

This  example  involves  data  on  the  survival  and  development  times  of  migratory  locusts  treated 
with  the  pathogens  Nosema  locustae  and  Nosema  cuneatum.  Although  these  disease  organisms  do  not 
inflict  as  much  mortality  as  would  be  desired  for  a  control  agent,  it  was  believed  that  they  might  provide 
additional  crop  protection  by  prolonging  development  through  reduced  consumption  and  activity.  The 
objective  of  this  study  was  to  determine  and  compare  the  effects  of  A'^.  locustae  and  A^.  cuneatum 
infections,  at  a  number  of  levels,  on  development  and  mortality  of  locusts. 

The  study  was  carried  out  using  individually  caged  third  instar  locusts  which  were  innoculated 
with  Nosema  spores  at  one  of  6  rates.  The  two  'failure'  types  were  death  and  the  attainment  of  the  adult 
stage.  The  time  to  reach  one  of  these  endpoints  was  recorded  for  each  locust.  The  study  was  stopped 
after  36  days,  and  any  locusts  still  alive  and  in  a  subadult  stage  were  considered  to  be  censored.  All  three 
competing  risks  models  were  applied  to  the  data. 
1.        The  mixture  model 

This  model  was  applied  to  the  locust  data  simply  by  computing  the  proportions  of 
individuals  dying  or  reaching  adulthood,  and  examing  the  observed  (conditional)  distributions  of 
death  and  development  times  (see  table  below).  This  is  appropriate  for  the  R  locustae  data 
because  there  was  virtually  no  censoring,  but  is  not  a  full  fitting  of  the  model  for  the  N.  cuneatum 
data.  This  approach  corresponds  to  the  customary  analysis  of  such  data. 

The  mean  observed  development  and  death  times  were  fairly  constant  over  all  dosage 
levels  of  spores  of  both  Nosema  species.  If  it  were  not  known  that  conditional  development  times 
can  be  misleading  about  the  effects  of  a  treatment  on  competing  failure  types,  one  would  likley 
conclude  that  Nosema  does  not  affect  development.  The  proportion  of  locusts  completing 
development  decreased  with  increasing  spore  dose,  but  these  proportions  do  not  tell  much  about 
the  impact  of  Nosema  on  development  because  such  a  decrease  could  simply  have  been  due  to 
the  increased  mortality  with  increasing  spore  doses. 
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Percent 


Ave.  Time  to 


Spore 

Species 

Dosage 

Num. 

Mort. 

Devel. 

Mort.  Devel. 

locustae 

0 

81 

4 

96 

26.3  27.5 

lO' 

82 

2 

98 

23.5  26.6 

lO' 

77 

39 

60 

29.0  28.4 

10" 

75 

72 

28 

26.8  26.6 

lO' 

89 

90 

9 

27.8  29.0 

10' 

88 

99 

1 

26.2  27.0 

cuneatum 

0 

75 

1 

99 

20.0  28.1 

lO' 

90 

3 

96 

19.0  28.7 

10" 

75 

3 

96 

26.5  27.5 

10" 

72 

32 

57 

28.1  29.3 

10^ 

87 

59 

11 

27.7  29.7 

10' 

86 

77 

5 

25.7  31.0 

2.        The  latent  failure  time  model 

As  can  be  seen  in  the  following  table,  the  estimated  means  of  the  distributions  of  latent 

development  and  death  times  showed  trends  with  increasing  spore  dose. 

Ave.  Latent 

Ave.  Latent 

Species 

Spore  Dosage 

Time  to  Death 

Time  to  Devel. 

locustae 

0 

40.5 

27.6 

10^ 

58.1 

26.7 

10^ 

31.9 

30.4 

10^ 

27.4 

29.5 

10^ 

28.2 

35.2 

26.2 

36.4 

cuneatum 

0 

117.1 

28.1 

10^ 

119.2 

28.7 

10^ 

48.1 

27.5 

10^ 

34.5 

31.7 

10* 

32.8 

42.0 

10' 

28.8 

41.9 
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If  the  assumption  of  independence  of  the  latent  times  were  true,  these  results  would 
indicate  that  as  spore  dose  increases,  expected  survival  time  of  the  locusts  decreases,  especially 
for  N.  locustae.  They  would  also  indicate  that  even  if  all  insects  were  to  survive,  they  would  take 
much  longer  to  complete  development  to  adulthood  under  the  influence  of  high  doses  (>  10*)  of 
Nosema  spores  of  either  species.  There  would  also  be  an  indication  that  A^.  cuneatum  prolongs 
development  to  a  greater  extent  than  N.  locustae. 
3.        The  cause-specific  hazard  model 

The  estimated  development  hazard  and  death  hazard  functions  demonstrated  conclusively 
that  both  tiie  mortality  and  development  processes  were  affected  by  Nosema.  The  deatii  hazard 
function  was  nonzero  for  A^.  locustae  doses  greater  than  or  equal  to  10^  and  for  A^.  cuneatum  doses 
greater  than  or  equal  10*.  Analysis  of  the  death  hazard  functions  via  the  proportional  hazards 
model  gave  the  following  proportionality  factors  for  doses  greater  tiian  or  equal  to  10*. 


Proportionality  Factor 

Spore  Dosage 

locustae  cuneatum 
10*                                      1.00  0.25 
10'                                     0.95  0.39 
10^                                     1.51  0.72 


The  development  hazard  function  seemed  unaffected  by  N.  locustae  at  a  dose  of  10^  spores,  but 
for  doses  of  10^  spores  or  more  tiie  hazard  function  decreased  toward  0  in  a  regular  manner.  For  A^. 
cuneatum,  the  development  hazard  function  was  unaffected  for  doses  less  than  or  equal  to  10^,  but  was 
just  as  depressed  as  for  A^.  locustae  at  higher  doses.  Proportional  hazards  analysis  of  the  development 
hazard  functions  (see  table  below)  showed  that  the  effects  of  equivalent  doses  of  A^.  locustae  and  A^. 
cuneatum  on  development  were  very  similar  except  at  the  lowest  doses,  10^  and  10^  spores. 


Proportionality  Factor 


Spore  Dosage 

locustae 

cuneatum 

0 

1.00 

0.89 

10^ 

1.40 

0.74 

10^ 

0.42 

1.02 

10* 

0.42 

0.32 

10' 

0.08 

0.06 

10^ 

0.02 

0.03 
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The  death  hazard  ftinctions  were  reanalysed  using  a  proportional  hazards  model  with  a 
time-dependent  covariate.  The  value  of  this  covariate  was  0  (zero)  if  the  locust  had  not  yet  achieved  the 
fifth  instar,  and  was  1  after  it  became  a  fifth-instar  locust.  Obviously,  this  covariate  contains  crude 
information  about  the  nearness  of  an  individual  to  the  adult  stage,  and  can  be  used  to  investigate  the 
relationship  between  the  two  competing  endpoints.  The  following  proportionality  factors  were  obtained: 


Proportionality  Factor 


Spore  Dosage 

Age 

locustae 

cuneatum 

10* 

inst<5 

1.00 

0.88 

inst=5 

0.79 

0.15 

10^ 

inst<5 

0.87 

0.95 

inst=5 

0.69 

0.16 

10^ 

inst<5 

1.51 

1.30 

inst=5 

1.19 

0.22 

This  analysis  shows  that  for  locusts  treated  with  A^.  cuneatum,  there  was  a  large  reduction  in  the 
death  hazard  upon  the  onset  of  the  fifth  instar,  whereas  for  those  treated  with  N.  locustae  the  reduction 
was  very  small.  Thus  for  A^.  cuneatum,  the  competing  endpoints  strongly  influence  each  other. 

This  analysis  also  suggests  that  the  modes  of  action  of  the  two  species  of  Nosema  are  quite 
different.  A  related  interesting  result  of  this  analysis  is  that  when  the  time-dependent  covariate  is  in  the 
model,  there  is  no  longer  a  significant  difference  between  the  species  of  Nosema.  Apparently  most  of  the 
difference  between  species  in  their  death  hazard  functions  can  be  explained  by  the  reduction  in  the  death 
rate  which  occurs  at  the  onset  of  the  fifth  instar  for  locusts  treated  with  N.  cuneatum.  The  use  of  the 
cause  specific  hazard  model  has  provided  information  that  no  other  procedure  would  have  detected,  and 
which  could  not  have  been  even  guessed  at  by  a  superficial  summary  of  the  data. 
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ABSTRACT 

The  statistical  procedure,  analysis  of  covariance  has  been  used  in  several  contexts.  The  most 
common  description  of  analysis  of  covariance  is  to  adjust  the  analysis  for  variables  that  could  not  be 
controlled  by  the  experimenter.  For  example,  a  researcher  can  remove  the  differential  effects  of  a  fertility 
trend  by  using  a  randomized  complete  block  design  structure,  but  it  may  not  be  possible  to  control  the 
number  of  plants  per  plot  of  land.  The  researcher  wishes  to  compare  varieties  as  if  each  plot  had  the  same 
number  of  plants.  The  analysis  of  covariance  is  a  procedure  which  can  compare  variety  means  by  first 
adjusting  for  the  differential  number  of  plants  per  plot. 

The  analysis  of  covariance  described  here  is  in  a  more  general  context  than  that  of  adjusting  for 
variation  due  to  uncontrollable  variables.  The  ANALYSIS  OF  COVARIANCE  is  defined  as  a  method 
for  comparing  several  regression  surfaces  or  lines,  one  for  each  treatment  or  treatment  combination,  where 
there  is  possible  a  different  regression  surface  for  each  treatment  or  treatment  combination. 

1.  Introduction 

The  experimental  situation  involves  randomly  assigning  experimental  units  to  treatment  i, 
applying  the  treatments  and  measuring  yi,Xiij,X2ij,...,Xqij  on  each  experimental  unit,  where 

yjj        is  the  dependent  measure, 

Xjij       is  the  first  independent  variable  or  covariate, 

X2ij       is  the  second  independent  variable  or  covariate, 

x^ij       is  the  kth  independent  variable  or  covariate. 
At  this  point,  the  experimental  design  is  a  one-way  treatment  structure  in  a  completely  randomized  design 
structure  with  k  covariates.  The  covariance  model  consisting  of  a  linear  function  of  the  covariates  or 
independent  variables  is 

Xij  =  Poi  +  PliXiij  +  p2iX2ij  +  ...  +  PkiXkij  +  Eij  (1.1) 

for  i=l,2,...t  where  t  is  the  number  of  treatments,  j  =  l,2,...,ni,  and  the  ~  N(0,  d^).  The  important  thing 
to  note  about  this  model  is  the  mean  of  the  y  values  for  given  values  of  the  x's  depends  on  the  values  of 
the  x's  as  well  as  on  the  treatment  or  treatment  combinafion  or  population  from  which  the  data  were 
collected.  The  analysis  of  covariance  is  a  strategy  for  making  decisions  about  the  form  of  the  covariance 
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model  through  testing  hypotheses  and  then  comparing  the  treatments  by  comparing  the  estimated  responses 
from  the  resulting  regression  models. 

Possible  forms  of  hypotheses  are: 

Hoi:  Phi  =  Ph2  =  •••  =  Pht  =  0  vs  H,i:  (not  H^i:),  that  is,  all  the  slopes  for  the  hth 
covariate  are  zero,  or 

H02:  Phi  =  Ph2  =  •••  =  Pht  vs  H^:  (not  H^jO,  that  is,  all  the  slopes  for  the  hth 
covariate  are  equal  meaning,the  surfaces  are  parallel  in  the  direction  of  the  hth  covariate. 

The  analysis  of  covariance  model  is  a  combination  of  the  analysis  of  variance  model  and  the 
regression  model.  An  experiment  is  designed  to  purchase  a  certain  number  of  degrees  of  freedom  for  error 
(generally  with  out  the  covariates)  and  the  experimenter  is  willing  to  sell  some  of  those  degrees  of 
freedom  for  good  covariates  which  will  help  reduce  the  magnitude  of  the  error  variance.  The  philosophy 
in  this  book  is  to  select  the  simplest  possible  expression  for  the  covariate  part  of  the  model  before  making 
treatment  comparisons.  This  process  of  model  building  to  determine  the  simplest  adequate  form  of  the 
model  follows  the  principle  of  parsimony  as  well  as  helps  guard  against  foolishly  selling  degrees  of 
freedom  for  error  to  obtain  unnecessary  covariate  terms  in  the  model.  Thus  the  strategy  for  analysis  of 
covariance  begins  with  testing  hypotheses  like  above  to  make  decisions  about  the  form  of  the  covariate 
or  regression  part  of  the  model. 

The  structure  of  the  following  chapters  leads  one  through  the  forest  of  analysis  of  covariance  by 
starting  with  the  simple  model  with  one  covariate  through  the  complex  process  involving  analysis  of 
covariance  in  split-plot  and  repeated  measures  designs.  Other  topics  discussed  are  multiple  covariates, 
experiments  involving  blocks,  and  graphical  methods  for  comparing  the  models  for  the  various  treatments. 

2.        One-way  Analysis  of  Covariance-One  Covariate  in  a  Completely  Randomized  Design 
Structure 

2.1      The  Model 

Suppose  there  are  N  homogeneous  experimental  units  which  are  randomly  divided  into  groups  of 
nj  units  per  group  where 

i  =N. 

i-l 

Each  of  the  t  treatments  of  a  one-way  treatment  structure  is  randomly  assigned  to  a  group  of  experimental 
units,  providing  a  one-way  treatment  structure  in  a  completely  randomized  design  structure.  Let  y^ 
(dependent  variable)  denote  the  jth  observation  from  the  ith  treatment  and  x^j  denote  the  covariate 
(independent  variable)  corresponding  to  the  (i,j)th  experimental  unit.  Assume  that  the  mean  of  y^  can  be 
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expressed  as  a  linear  function  of  the  covariate  with  possibly  different  linear  functions  required  for  each 
treatment.  It  is  important  to  note  that  the  mean  of  an  observation  from  the  ith  treatment  group,  depends 
on  the  particular  treatment  as  well  as  on  the  values  of  the  covariate  (independent  variable). 

The  analysis  of  covariance  model  for  a  one-way  treatment  structure  with  one  covariate  in  a 
completely  randomized  design  structure  is 

yy  =  oti  +  piXij  +  Eij, 
i  =  l,2,...,t.  j  =  l,2,...,n  (2.1) 
where  the  mean  of  y■^  for  a  given  value  of  X  is     |x  =  otj  +  Pi^  and  for  making  inferences,  it  is  assumed 
that 

ejj  ~  iid  N(0,cr^). 

Model  (2.1)  has  t  intercepts,  ai,...,ai,  t  slopes  Pi,...,pt  and  one  variance  o^.  Before  one  can  continue  with 
the  analysis  of  this  model,  one  must  be  sure  that  the  data  from  each  treatment  can  in  fact  be  described  by 
a  simple  linear  regression  model.  Various  regression  diagnostics  should  be  run  on  the  data  before 
continuing.  The  equal  variance  assumption  should  also  be  checked.  If  the  simple  linear  regression  model 
is  not  adequate  to  describe  the  data  for  each  treatment  another  model  must  be  selected  before  continuing 
with  the  analysis  of  covariance.  The  analysis  of  covariance  involves  comparing  the  t  slopes  and 
comparing  the  distances  between  the  regression  lines  (surfaces)  at  preselected  values  of  X.  The  analysis 
of  covariance  computations  are  typically  presented  in  summation  notation  with  little  emphasis  on 
interpretations.  In  this  paper  the  various  covariance  models  are  expressed  in  terms  of  matrices  (see 
Chapter  6  of  Milliken  and  Johnson  (1984))  and  their  interpretations  are  discussed.  The  computer  is  used 
as  the  mode  of  doing  the  analysis  of  covariance  computations. 
The  matrix  form  of  Model  (2.1)  is 


>'ll 

1 

^11 

0 

0  • 

.•  0 

0~ 

a, 

1 

In, 

0 

0  . 

0 

0 

Pi 

3^21 

0 

0 

1 

^21  • 

•  0 

0 

0 

0 

1 

•  0 

0 

P2 

y.j 

0 

0 

0 

0  . 

.  1 

a, 

0 

0 

0 

0  • 

•  1 

P,. 
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which  is  expressed  in  the  form  of  a  Hnear  model  as  ^  =  2£P  +  £•  The  vector  ^  denotes  the  collection  of 
slopes  and  intercepts,  the  matrix  X  is  the  design  matrix  and  the  vector  e  represents  the  random  errors. 

2.2  Estimation 

The  least  squares  estimator  of  the  parameter  vector  ^  is 

^  =  (x'x)-^xV. 

But  the  least  squares  estimator  of  p  can  also  be  found  by  obtaining  the  least  squares  estimator  of  each  pair 
of  parameters  (Oj,  ft)  by  fitting  the  simple  linear  model  to  the  data  from  each  treatment.  For  data  from 
the  ith  treatment,  fit  the  model 


1\ 

}  \ 

k 

which  is  expressed  as  ^  =  X3;  +  e^.  The  least  squares  estimator  of  ^  is 

the  same  as  the  estimator  obtained  for  a  simple  linear  regression  model.  The  estimates  of  ft  and  (x^  in 
summation  notaUon  are 

n. 

E  x.  y  .-n.x.y. 
0,  =  £L  

E  Xij  -n^x, 
j-i 

and 

&,  =  y,  -|5,x,. 
The  residual  sum  of  squares  for  the  ith  model  is 

SSRes,  =  E  (y,.  -  a,  - 

which  is  based  on  nj-2  degrees  of  freedom.  After  testing  the  equality  of  the  treatment  variances,  the 
residual  sum  of  squares  for  model  (2.1)  can  be  obtained  by  pooling  residual  sums  of  the  squares  for  each 
of  the  t  models,  i.e.,  sum  them  together  to  obtain 

SSRes  =  E  SSReSj.  (2-4) 
The  pooled  residual  sum  of  squares,  SSRes,  is  based  on  the  pooled  degrees  of  freedom 
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±f.,,^^  =  Ein,-2)-En-2t  =  N  -2t. 


The  best  estimate  of  the  experimental  unit  variance  is 

6^  =  SSRes/(N  -  2t). 

The  sampling  distribution  of 

(N  -  2t)dVo^ 

is  central  chi-square  with  (N  -  2t)  degrees  of  freedom.  The  sampling  distribution 
estimator, 

is  normal  with  mean  ^  =  (ai,pi,...,at,pj)  and  variance-covariance  matrix  cf-(XX)'\  which  can  be  written 


as 


(2c;x)-^  0 
0      ..•  (X'X)-^ 


(2.5) 


where 


1  X. 


1  X 


2.3      Strategy  for  Determining  the  Form  of  the  Model 

The  main  objective  of  an  analysis  of  covariance  is  to  compare  the  t  regression  lines  at  several 
predetermined  fixed  values  of  the  covariate,  X.  Depending  on  the  values  of  the  pj,  there  are  various 
strategies  that  one  takes  when  comparing  the  regression  lines. 

The  first  question  which  needs  to  be  answered  is,  does  the  mean  of  y  given  X  depend  on  the  value 
of  X?  That  question  can  be  answered  "statistically"  by  testing  the  hypothesis  Hq^:  E(y-- 1  X  =  x)  =  ttj  vs. 
H31:  E(yij  I  Xjj  =  Xjj)  =  ot,  +  pjX  for  i  =  l,2,...,t.  The  hypothesis  is  equivalent  to 

Hoii  pj  =  p,  =  ...  =  p^  =  0  vs.  H,i:  (not  Ho).  (2.6) 
The  null  hypothesis  states  that  none  of  the  treatments  have  means  which  depend  linearly  on  the  value  of 
the  covariate,  X. 

The  principle  of  conditional  error  (Milliken  and  Johnson,  (1984))  or  model  comparison  method 
(Draper  and  Smith  (1981))  provides  an  excellent  way  of  obtaining  the  desired  test  statistic.  The  model 
restricted  by  the  conditions  of  the  null  hypothesis,  Hqi,  is 
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y.,  =  +  e..  i  =  i,2,...,t,  j  =  1,2,... A-  (2.7) 
Model  (2.7)  is  the  usual  analysis  of  variance  model  for  the  one-way  treatment  structure  in  a  completely 
randomized  design  structure.  The  residual  sum  of  squares  for  model  (2.7)  is 

SSRes(Hoi)  =  E  E(y..  -  y,}^  (2.8) 

i-l  j-1 


which  is  based  on 


d.f.  ssRcs(H„)  =  N  - 1  degrees  of  freedom 


(where  the  mean  of  the  model  under  Hqi  has  t  parameters).  The  sum  of  squares  due  to  deviations  from 
Hqi,  denoted  by  SSHqi  is, 

SSHoi  =  SSRes(Hoi)  -  SSRes,  (2.9) 

which  is  based  on 

^•^•ssRes(H„)  "  ^'f-  SSRes  =  (N  -  t)  -  (N  -  2t)  =  t  dcgrccs  of  freedom. 

The  sampling  distribution  of  SSHqi/o^  is  a  noncentral  chi-square  distribution  with  t  degrees  of  freedom 
where  the  non-centrality  parameter  is  zero  if  and  only  if  Hqi  is  true.  A  statistic  for  testing  Hqi  versus  H^i 
is 

F     =  ^^"q^^^  (2.10) 

Ho,  ^ 

and,  when  H^i  is  true,  the  sampling  distribution  of  F„  is  a  central  F  distribution  with  t  and  N  -  2t  degrees 
of  freedom. 

If  one  fails  to  reject  Hqi,  then  one  can  conclude  that  the  means  of  the  treatments  do  not  depend 
linearly  on  the  value  of  the  covariate,  X.  In  this  case,  the  next  step  in  the  analysis  is  to  use  analysis  of 
variance  to  make  comparisons  between  treatment  means,  i.e.,  compare  the  0^,  i  =  l,2,...,t  (see  chapter  1, 
Milliken  and  Johnson  (1984)).  Recall  it  has  already  been  determined  that  the  simple  linear  regression 
model  adequately  describes  the  data.  Thus  if  the  slopes  are  zero,  then  conclude  the  models  are  of  the  form 
Yij  =  oq  +  Eij. 

If  Hqi  is  rejected,  then  conclude  that  the  mean  of  y  does  depend  linearly  on  the  value  of  the 
covariate  X  for  at  least  one  of  the  treatments.  In  this  case,  the  next  step  in  the  analysis  of  covariance  is 
to  determine  whether  or  not  the  means  of  the  treatments  depend  on  the  covariate  X  differently  (as 
represented  by  unequal  slopes  which  provide  non-parallel  lines).  A  test  for  homogeneity  or  equality  of 
the  slopes  will  answer  that  question.  The  appropriate  null  hypothesis  is 

H02:  E(yij  I  X  =  x)  =  ttj  +  (3x  vs.  H^:  Eiy-^^  |  X  =  x)  =     +  p^x  or  equivalently 

Ho^:  Pi  =      =  ...  =     =  [3  vs.        (not  H^^)  (2.11) 
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where  p  is  unspecified  and  represents  the  common  slope  of  the  t  parallel  regression  lines.  The  model  in 
the  form  of  (2.1)  which  satisfies  the  conditions  of  is 

=  +  pXij  +  i  =  l,2,...,t,  j  =  l,2,...,n,  (2.12) 
which  represents  t  parallel  lines  each  with  slope  p  and  intercepts  ai,...,at.  The  matrix  form  of  model 
(2.12)  is 

■"l  0  ...  0  X, 


J  in, 


1  0      0  X 


In, 


0  1  ...  0  X21 

0  1  ...  0  x^ 

0  0  ...  1  x^j 

0  0  ...  1  x„ 


a. 


a. 


(2.13) 


The  residual  sum  of  squares  for  model  (2.13)  is 


SSRes  (H02)  =  E  E  (y,  -  a,  -  ^xS 


(2.14) 


'I  j-i 


where 


a.  i  =  1,2  t  and  ^ 

denote  the  least  squares  estimators  of  the  corresponding  parameters  from  model  (2.13).  The  residual  sum 
of  squares  in  (2.14)  is  based  on 

d.f.  ssRcs(Hj  ~        degrees  of  freedom 
as  the  mean  of  model  (2.12)  has  t  +  1  parameters.  The  sum  of  squares  due  to  deviations  from  H02  is 

SSH02  =  SSRes(Ho2)  -  SSRes  (2.15) 

which  is  based  on 

d-f-  ssRes(H„)  -  d.f.  ssRcsiuj  =  ^ "  1  ^^egrecs  of  freedom. 
The  sampling  distribution  of  SSH^/d^  is  non-central  chi-square  with  t  -  1  degrees  of  freedom  where  the 
non-centrality  parameter  is  zero  if  and  only  if      is  true.  The  statisfic  used  to  test  is 

SSHJt  -  1 


6^ 


(2.16) 
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which  has  sampling  distribution  F  with  t  -  1  and  N  -  2t  degrees  of  freedom.  If  one  fails  to  reject  H^,  then 
conclude  that  the  lines  are  parallel  (equal  slopes)  and  proceed  to  compare  the  distances  between  the 
parallel  regression  lines  by  comparing  their  intercepts,  a^'s  (which  is  discussed  in  Section  2.4).  Figure  2.1 
displays  the  relationships  between  the  treatment  means  as  a  function  of  the  covariate  X  when  the  lines  are 
parallel.  Since  the  lines  are  parallel,  i.e.,  the  distance  between  any  two  lines  is  the  same  for  all  values  of 
X,  a  comparison  of  the  intercepts  is  a  comparison  of  the  distances  between  the  lines. 
If  one  rejects  H02,  then  conclude  that  at  least  two  of  the  regression  lines  have  unequal  slopes  and  hence, 
the  set  of  lines  are  not  parallel.  Figure  2.2  displays  a  possible  relationship  between  the  means  of 
treatments  as  a  linear  function  of  the  covariate  for  the  non-parallel  line  case.  When  the  lines  are  not 
parallel,  the  distance  between  two  lines  depends  on  the  value  of  X,  thus  nonparallel  line  case  is  called 
covariate  by  treatment  interaction. 

2.4      Comparing  the  Treatments  or  the  Regression  Lines 

The  correct  method  for  comparing  the  distances  between  the  regression  lines  depends  on  the 
decision  one  made  concerning  the  slopes  of  the  models.  If  the  experimenter  rejects  Hqi  in  (2.6)  and  fails 
to  reject  in  (2.11),  the  resulting  model  is  a  set  of  parallel  hnes  (equal  slopes).  A  property  of  two 
parallel  lines  is  that  they  are  the  same  distance  apart  for  every  value  of  X.  Thus,  the  distance  between 
any  two  lines  can  be  measured  by  comparing  the  intercepts  of  the  two  lines.  When  the  lines  are  parallel, 
contrasts  between  the  intercepts  are  used  to  compare  the  treatments.  When  the  slopes  are  unequal  there 
are  two  types  of  comparisons  that  are  of  interest  to  the  experimenters,  namely,  comparing  the  distances 
between  the  various  regression  lines  at  several  values  of  the  covariate  X  and  comparing  specific 
parameters,  such  as  slopes  or  the  models  evaluated  at  the  same  or  different  values  of  X. 
2.4.1    Equal  Slope  Model 

At  this  step  in  the  analysis  remember  that      was  not  rejected,  thus  the  model  used  to  describe 
the  mean  of  y  as  a  function  of  the  covariate  is  of  the  form 

y^- a,  +  pX,j  =  e^j.  (2.17) 
The  residual  sum  of  squares  for  model  (2.17)  is  SSRes(Ho2)  which  was  given  in  equation  (2.13).  The  first 
hypothesis  to  be  tested  is  that  the  distances  between  the  lines  are  equal,  which  is  equivalent  to  testing  the 
hypothesis  that  the  intercepts  are  equal  as 

H03:  tti  =  02  =  ...  =  ot^  =  a  vs.  H^:  (not  H^)  (2.18) 
where  a  is  unspecified.  The  model  in  (2.17)  restricted  by  the  conditions  of  H^^  is 
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Plot  of  Y*X.     Symbol  used  is  ' +'  . 
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Figure  2. 1    Graph  of  parallel  lines  models— common  slopes. 
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Plot  of  Y*X.     Symbol  used  is  ' +'  . 
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Figure  2.2  Graph  of  non-parallel  lines  models — unequal  slopes. 
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y,^  =  a+pX,^  +  e,^.  (2.19) 
Model  (2.19)  is  a  single  simple  linear  regression  model  which  is  to  be  fit  to  aU  of  the  data  and  the 
corresponding  residual  sum  of  squares  is 

SSRes(Ho3)  =  i:  E  (y..  -  d  -      )^  (2.20) 

i-l  j-1 

which  is  based  on 

d.f.  ssRcs(H,^  =  N  -  2  degrees  of  freedom 

where  &  and  3  are  least  squares  estimators  of  a  and  p  from  model  (2.19). 

The  sum  of  squares  due  to  deviation  from  Hq,,  given  that  the  slopes  are  equal  (Hq^  is  true),  is 
SSHo3  =  SSRes(Ho3)  -  SSRes(Ho2)  (2.21) 

which  is  based  on 

d-f-  ssRescHj  -  d.f.  ssRes(Hj)  =  ^ "  ^  ^cgrecs  of  freedom. 


The  appropriate  test  statistic  is 

^"   ~  SSRes(H^)/(N  -  t  -  1) 


The  sampling  distribution  of  F^^  is  that  of  a  non-central  F  distribution  with  t  -  1  and  N  -  t  -  1 
degrees  of  freedom.  If  H03  is  not  rejected  conclude  that  all  of  the  data  comes  from  a  single  regression 
model  with  slope  p  and  intercept  a,  i.e.,  there  are  no  treatment  differences.  If  H03  is  rejected,  then 

conclude  that  the  distances  between  one  or  more  pairs  of  lines  are  different  from  zero.  Since  the  lines  are 
parallel,  the  distance  between  any  two  lines  can  be  compared  at  any  chosen  value  of  X.  If  the  distance 
between  any  two  lines  is  compared  at  X  =  0,  it  is  a  comparison  between  the  difference  of  two  intercepts 
as  a  -  a..  A  multiple  comparison  procedure  can  be  used  to  compare  the  distances  with  an  LSD  type 
being  used  for  controlling  the  comparison  wise  error  rate  and  a  Bonferroni  or  Scheffe  type  being  used  to 
control  the  experiment  wise  error  rate.  To  use  a  multiple  comparison  procedure,  one  must  compute  the 
standard  error  of  the  difference  between  two  a's  as 

S,^,  =  [S,^  .  S,^,  -  2Cov(&„  <3^f  (2.23) 

where 

Sft  is  the  variance  of  6l.  and  Cov(d,  Cl.,)  is  the  covariance  between  ft.  and 

If  there  are  specific  planned  comparisons  (such  as  linear  or  quadrafic  effects  for  treatments  with 
quanfitative  levels)  between  the  treatments,  those  comparisons  would  be  made  by  construcdng  the 
necessary  contrasts  between  the  intercepts. 
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When  analysis  of  covariance  was  first  developed,  it  was  mainly  used  to  adjust  the  mean  of  y  for 
a  selected  value  of  the  covariate.  The  value  usually  selected  was  the  mean  of  the  covariate  from  all  t 
treatments.  Thus  the  term  adjusted  means  was  defined  as  the  mean  of  y  evaluated  at 

X  =  X...,  where  x.. 

is  the  mean  value  of  all  the  x^j's.  The  estimators  of  the  means  of  treatments  at 

X=x.., 

called  the  adjusted  means,  are 

|1y,Ipx.p;  =  6.,  *Px..  i  =  l,2....,t.,  (2.24) 

where  Cl.  and  (5  are  least  squares  estimators  of  oq  and  P  from  model  (2.18).  The  covariance  matrix  of 
the  adjusted  means  can  be  constructed  from  the  elements  of  the  covariance  matrix  of  6l^,....6l^  and  3- 

The  standard  errors  of  the  adjusted  means  are  computed  as 


Sl  +  X^S.^  +  2Xcov((l,P) 


,  i  =l,2,...,t. 


One  hypothesis  of  interest  is  that  of  the  equality  of  the  adjusted  means.  This  hypothesis  can  be  expressed 
as 

^04^   Mvjpx-pi  =  -  =  Mvjpx.pi  vs.   H(not  HJ. 
But  since  the  lines  are  parallel,  the  difference  between  two  adjusted  means  is  the  difference  between 
intercepts  as 

Pvjpx  »  PJ  "  ^^YjpX  -  pi  ^  ^1  " 

thus      is  equivalent  to  H03. 

Preplanned  treatment  comparisons  and  multiple  comparison  procedure  can  be  carried  out  to 
compare  the  adjusted  means  by  computing  the  standard  error  of  the  difference  between  pairs  of  adjusted 
means.  Since  the  difference  of  two  such  means  is  -  the  standard  error  in  (2.23)  can  be  used  with 
any  selected  multiple  comparison  procedure.  Contrasts  between  adjusted  means,  which  are  also  contrasts 
between  the  intercepts,  measuring  linear,  quadratic,  etc  trends  should  be  used  when  appropriate. 
2.4.2    Unequal  Slope  Model-Covariate  bv  Treatment  Interaction 

When  H02  is  rejected,  it  is  concluded  that  the  non-parallel  lines  model  (2.1)  is  necessary  to 
adequately  describe  the  data.  The  graph  of  such  a  possible  situation  was  given  in  figure  (2.2).  Since  the 
lines  are  not  parallel,  the  distance  between  any  two  lines  depends  on  which  value  of  the  covariate  is 
selected.  This  is  called  covariate  by  treatment  interaction.  In  the  nonparallel  lines  case  a  comparison  of 
the  intercepts  is  only  a  comparison  of  the  lines  at  X  =  0.  That  will  generally  be  a  meaningful  comparison 
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only  when  X  =  0  is  included  in  or  is  close  to  the  range  of  X  values  in  the  experiment.  The  equal  intercept 
hypothesis  given  the  slopes  are  unequal  is  expressed  as 

H05:  E(yij  I  X  =  xjj)  =  a  +  pjXy  vs  H^:  E(yij  |  X  =  x^)  =     +  piXy  or  equivalenUy, 
H05:  ttj  =  02  =  ....  =  Ot  =  a  given  (3;    pj-  vs  H„(not  Hq). 
The  model  comparison  method  is  used  to  test  this  hypothesis.  Model  (2.1)  restricted  by  the  conditions 
of  is 

y,  =  a+ftX,  +  e,  (2.25) 
and  the  corresponding  residual  sum  of  squares  as 

SSRes(Ho3)  =  i:  i:(y,.  -  &  -  3.X..)2  (2-26) 

i-l  j-1 

which  is  based  on 

d.f.  ssRes(H,^  =  N  - 1  -  1  degrees  of  freedom. 

The  values  of  A  and  Pj  in  (2.26)  are  the  least  squares  estimators  of  the  parameters  of  model  (2.25).  The 

sum  of  squares  due  to  deviations  from  H05  is 

SSH05  =  SSRes(Ho5)  -  SSRes, 

which  is  based  on 

^•f-  ssRes(Hj  -  d-f-  SSRes  =  t "  1  dcgrccs  of  frecdom. 

The  test  statistic  is 

F     =    SSH,3/(t  -  1)  ^2.27) 
SSRes/(n  -  2t)' 


The  conclusion  one  makes  at  X  =  0  may  be  different  from  that  made  at 

X  =  X..  or  X  =  Xo, 

where  Xq  is  some  other  fixed  value  of  the  covariate.  Suppose  the  experimenter  wants  to  compare  the 
distances  between  the  lines  at  a  selected  value  of  X,  say  X  =  Xq.  The  hypothesis  to  be  tested  is 

Py,Ip,x-p,x,  ^  MyJpx-p.x,  ^  •••  Py,|p.x-p.x,  ^  ^^x„  '  (2.28) 

where 

The  model  in  2. 1  can  be  equivalently  expressed  as 

y„  =  a  *  |i,X„  -  p,X„  .  p,X,,  .  e,.  =        .  ^^^^^ 
-  Pi(X„  -       *  e,. 
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The  model  restricted  by  is 

y,,  =Mvip,x.j,x.  *Pi(x,  -x^  -e,;  (2.30) 

and  the  corresponding  residual  sum  of  squares  is 

SSRes(HJ  =  E  E(y,.  -  p    -  ^.(X..  -  X^f,  (2-31) 

i-1  j-1 

which  is  based  on 

d-f-  ssRes(Hj  -  d-f-  ssRcs  =  N  - 1  -  1  dcgrccs  of  freedom. 
The  sum  of  squares  due  to  deviations  from  H^^  is  SSHq^  =  SSResCHo^)  -  SSRes,  which  is  based  on 


d  f-  ssRes(Hj  -  d-f-  ssR«  =  ^ "  ^  ^^^grees  of  freedom. 


The  resulting  test  statistic  is 


F    =     ^^^06^^     ^  (2.32) 
SSRes/(N  -  2t)' 


It  is  important  for  the  experimenter  to  make  comparisons  between  the  lines  at  several  different 
values  of  the  covariate.  The  usual  comparison  of  adjusted  means,  i.e.,  at 

X  =  X.. 

is  only  one  of  many  comparisons  which  are  probably  of  interest.  Figure  2.3  shows  three  possible  values 
of  X  (covariate)  at  which  one  might  make  comparisons.  The  corresponding  test  statistics  can  be  obtained 
by  expressing 

with  Xq  =  Xj,  Xq  =  X..  and  Xq  =  X* 
respectively.    If  a  hypothesis  corresponding  to  a  selected  value  of  Xq  is  rejected,  then  a  multiple 
comparison  procedure  could  be  used  to  determine  which  distances  between  pairs  of  lines  are  not  equal  to 
zero  or  the  preplanned  treatment  comparisons  could  be  made  at  Xq.  The  standard  error  of  the  difference 
between  two  lines  at 

is 

(2.33) 


1^ 


where  the  standard  errors  S.  can  be  obtained  from  the  standard  errors  of  the  intercept  parameters 
in  model  (2.29).  These  standard  errors  are  computed  with  the  assumption  the  two  models  have  no 
common  parameters  so  the  covariance  between  the  two  adjusted  means  is  zero.   Again,  preplanned 
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comparisons  and  LSD,  Bonferroni  or  Scheffe  types  of  multiple  comparison  procedures  can  be  used  to  help 
interpret  the  results  of  the  analysis  of  covariance.  For  most  experiments,  comparisons  should  be  made  for 
at  least  three  values  of  X,  one  in  the  lower  range,  one  in  the  middle  range  and  one  in  the  upper  ranges 
of  the  X's  obtained  for  the  experiment. 

The  experimenter  is  often  interested  in  determining  which  treatment  mean  responds  the  most  to 
a  change  in  the  value  of  the  covariate.  In  this  case  an  LSD  approach  can  be  used  to  make  size  a 
comparison  wise  tests  about  the  ft's.  Alternatively,  one  could  also  use  a  Bonferroni  or  Scheffe-type 
approach  to  control  the  experiment  wise  error  rate.  In  any  case,  the  standard  error  of  the  difference 
between  two  slopes  is 

where  Sj^  denotes  the  standard  error  associated  with  and  it  is  assumed  the  covariance  between  the  two 
parameters  is  zero  (which  is  the  case  here  since  the  two  models  do  not  have  common  parameters.)  The 

degrees  of  freedom  for  the  appropriate  percentage  point  is  N  -  2t.  Preplanned  treatment  comparisons  can 

be  made  by  comparing  the  slopes  of  the  various  models.   Section  2.7  shows  how  to  carry  out  the 

computations  via  SAS. 


2.5      Confidence  Bands  About  the  Difference  of  Two  Treatments 

When  the  slopes  are  unequal,  it  is  often  useftil  to  determine  the  region  of  the  covariate  where  the 
two  treatments  produce  significantly  different  responses.  A  confidence  band  can  be  constructed  about  the 
difference  of  two  treatment  models  and  the  region  of  the  covariate  where  the  confidence  band  does  not 
contain  zero  is  determined.  A  Scheffe  type  confidence  statement  should  be  used  to  provide 
experimentwise  error  protection.  The  difference  between  two  lines  for  treatment  1  and  2  at  X  =  Xq  is 

which  has  standard  error 
where 

An  example  of  the  construction  and  use  of  the  confidence  band  about  the  difference  of  two  regression 
lines  is  in  Section  2.8. 
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Plot  of  Y*X.     Symbol  used  is 
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Figure  2.3       For  non-parallel  lines  models  the  comparisons  between  the  treatments  depends  of  the 
selected  value  of  X. 
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2.6      Summary  of  Strategies 

Sections  2.3  and  2.4  describe  the  strategies  for  determining  the  form  of  the  model  and  the  resulting 
analyses.  Table  2.6. 1  lists  the  paths  for  model  determination,  but  there  are  possible  exceptions. 


Table  2.6.1      Strategy  for  determining  the  form  of  the  analysis  of  covariance  model  involying 
one  covariate  assuming  a  simple  linear  regression  model  will  describe  each 
treatment  data. 


a)  Test  the  hypothesis  that  the  slopes  are  zero 

i)  If  fail  to  reject,  compare  the  treatments  via  analysis  of  variance, 

ii)  if  reject  go  to  b) 

b)  Test  the  hypothesis  that  the  slopes  are  equal 

i)  if  fail  to  reject,  use  a  parallel  lines  model  and  compare  the 
treatments  by  comparing  the  intercepts  or  adjusted  means 
(LSMEANS). 

ii)  if  reject  to  c) 

c)  Use  the  unequal  slope  model  and 

i)  compare  the  slopes  of  the  treatments  to  see  if  treatments  can  be 
grouped  into  groups  with  equal  slopes 

ii)  compare  the  models  of  at  least  three  values  of  the  covariate,  low, 
middle  and  high  value 

iii)  construct  confidence  bands  about  the  difference  of  selected  pairs  of 
models. 


There  are  at  least  two  possible  exceptions  to  the  strategy.  First,  it  is  possible  to  reject  the  equal  slope 
hypothesis  and  fail  to  reject  the  zero  slope  hypothesis  when  there  are  both  positive  and  negative  slopes. 
In  this  case,  use  the  non-parallel  lines  model.  Second,  it  is  possible  to  fail  to  reject  the  zero  slopes 
hypothesis  when  in  fact  the  slope  of  the  parallel  lines  model  is  significantly  different  from  zero.  Here 
there  is  not  enough  information  from  the  individual  treatments  to  say  the  slopes  are  different  from  zero, 
but  the  combined  information  does  detect  the  linear  relationship.  In  this  case  use  the  common  slope  or 
parallel  lines  model. 

2.7      Analysis  of  Covariance  Computations  via  the  SAS  System 

The  SAS  System  can  be  used  to  compute  the  various  estimators  and  tests  of  hypotheses  discussed 
in  the  previous  sections.  The  SAS  System  statements  required  for  each  part  of  the  analysis  are  presented 
in  this  section  and  detailed  examples  are  discussed  in  Sections  2.8  and  2.9. 

AU  the  following  models  will  be  fit  assuming  that  the  data  were  read  in  by  the  following  statements: 
DATA  ONECOV;INPUT  TRT  Y  X; 
CARDS; 

The  required  SAS  System  statements  needed  to  fit  model  (2.1)  are 
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PROC  GLM;  CLASSES  TRT; 

MODEL  Y  =  TRT  X*TRT  /NOINT  SOLUTION; 
The  term  TRT  with  the  no-intercept  (NOINT)  option  generates  the  part  of  the  design  matrix 
corresponding  to  the  intercepts  and  enables  one  to  obtain  the  estimators  of  the  intercepts.  The  term 
X*TRT  generates  the  part  of  the  design  matrix  corresponding  to  the  slopes.  The  SOLUTION  option  is 
used  so  that  the  estimators  and  their  standard  errors  are  printed.  (PROC  GLM  does  not  provide  the 
estimators  when  there  is  a  CLASS  variable  unless  SOLUTION  is  specified.)   The  sum  of  squares 

corresponding  to  ERROR  is  SSRes  of  equation  (2.4)  and  the  MEAN  SQUARE  ERROR  is  ^^  the  estimate 
of  the  sampling  variance. 

The  type  III  (or  Type  IV)  sum  of  squares  corresponding  to  X*TRT  tests  Hqi  of  (2.6). 
The  Type  III  (or  Type  IV)  sum  of  squares  corresponding  to  TRT  tests  H05,  i.e.,  ai  =  ot2=...  =  (Xt  =  0 
given  that  the  unequal  slopes  are  in  the  model.  This  hypothesis  is  often  not  of  interest,  but  a  test  is 
available  in  case  the  experimenter  has  a  situation  where  a  zero  intercept  hypothesis  is  interpretable. 

Next,  to  test  H02  of  (2.11),  the  required  SAS  System  statements  are: 
PROC  GLM;  CLASSES  TRT; 
MODEL  Y  =  TRT  X  X*TRT  /  SOLUTION; 
The  type  III  sum  of  squares  corresponding  to  X*TRT  tests  H^.    The  type  III  sum  of  squares 
corresponding  to  X  tests  if  the  average  value  of  the  slopes  is  zero  and  the  type  III  sums  of  squares 
corresponding  to  TRT  tests  H05.   By  including  X  and/or  removing  the  NOINT  option,  the  model  is 
singular  and  the  provided  least  squares  solutions  is  not  directly  interpretable.  The  obtained  least  squares 
solution  satisfies  the  set-to-zero  restrictions  (see  chapter  6  of  Milliken  and  Johnson  (1984)).  If  one  uses 
the  model  statement  Y  =  TRT  X*TRT  X,  where  X*TRT  is  listed  before  X,  the  type  I  sum  of  squares 
corresponding  to  X*TRT  tests  Hqi  while  the  type  III  sum  of  squares  tests  H02.  A  list  of  type  I  and  III 
estimable  functions  can  be  obtained  and  used  to  verify  the  hypothesis  tested  by  each  sum  of  squares. 
If  one  fails  to  reject  H02,  the  parallel  lines  or  equal  slope  model  of  2.12  should  be  fit  to  the  data.  The 
appropriate  SAS  System  statements  are: 

PROC  GLM;  CLASS  TRT; 

MODEL  Y  =  TRT  X  /SOLUTION;. 

The  type  III  sum  of  squares  corresponding  to  TRT  is  SSH03  of  (2.21),  and  the  resulting  F  ratio 
tests  that  the  distances  between  the  lines  are  zero  given  that  the  parallel  line  model  is  adequate  to  describe 
the  data. 

Estimates  of  the  mean  of  y  given  X  =  x..  for  each  treatment,  which  are  often  called  adjusted 
means,  can  be  obtained  by  including  the  statement 
LSMEANS  TRT  /  STDERR  PDIFF; 
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after  the  MODEL  statement.  This  provides  the  adjusted  means 


the  estimate  of  the  treatment  mean  at  X  =  x  or  least  squares  mean,  for  each  treatment  and  the  option 
STDERR  provides  the  corresponding  standard  errors  of  the  adjusted  means.  The  PDIFF  option  provides 
significance  levels  for  t-tests  of 


(The  TDIFF  option  provides  the  values  of  the  t-tests  from  which  the  PDIFF  values  are  obtained.)  for  each 
pair  of  adjusted  means.  A  comparison  of  adjusted  means  is  also  a  comparison  of  the  o^j's  for  the  parallel 
lines  model.  The  significance  probabilities  can  be  used  to  construct  a  LSD  or  a  Bonferroni  multiple 
comparison  procedure  for  comparing  the  distances  between  pairs  of  lines. 

Any  comparisons  between  parameters  can  be  made  by  using  the  ESTIMATE  or  CONTRAST 
statement  in  the  GLM  procedure.  There  are  two  situations  where  such  statements  are  needed. 
First,  if  the  conclusion  is  that  the  slopes  are  not  equal,  then  one  can  apply  a  multiple  comparison 
procedure  in  order  to  compare  some  or  all  pairs  of  slopes.  This  is  easily  done  by  including  an 
ESTIMATE  statement  following  the  MODEL  statement  for  each  comparison  of  interest.  For  example, 
if  there  are  three  treatments  and  it  is  of  interest  to  compare  all  pairs  of  slopes,  then  the  following 
statements  would  be  used: 

PROC  GLM;  CLASSES  TRT; 

MODEL  Y=TRT  X*TRT  /  NOINT  SOLUTION; 

ESTIMATE  'B1-B2'  X*TRT  1  -1  0; 

ESTIMATE  'B1-B3'  X*TRT  1    0  -1; 

ESTIMATE  'B2-B3'  X*TRT  0    1  -1; 
Each  ESTIMATE  statement  produces  an  estimate  of  the  linear  combination  of  parameters,  a  computed  t- 
value,  and  it's  significance  level  using  the  residual  mean  square  as  the  estimate  of  cr^.  The  significance 
levels  obtained  from  these  comparisons  can  be  used  to  construct  a  LSD  or  Bonferroni  multiple  comparison 
procedure. 

For  the  unequal  slope  model,  adjusted  means  needs  to  be  obtained  X  =  x..  and  at  other  specified 
values  of  X.  Such  adjusted  means  can  be  used  to  make  comparisons  between  the  treatments.  The  SAS 
system  code  for  using  ESTIMATE  statements  to  obtain  adjusted  means  at  X  =  7.3  is 

ESTIMATE  Tl  AT  7.3'  TRT  1  0  0  X*TRT  7.3  0  0; 

ESTIMATE  'T3  AT  7.3'  TRT  0  1  0  X*TRT  0  7.3  0; 

ESTIMATE  'T3  AT  7.3'  TRT  0  0  1  X*TRT  0   0  7.3;. 
Contrasts  between  the  adjusted  means  at  X  =  7.3  can  be  obtained  by  subtracting  the  respective  estimate 
statement  values  from  the  adjusted  mean  estimate  statements  as 
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ESTIMATE  'Tl  -  T2  AT  7.3'  TRT  1  -1  -  X*TRT  7.3  -7.3  0; 

ESTIMATE  'Tl  -  2T  +  T3  AT  7.3'  TRT  1  -2  1  X*TRT  7.3  -14.6  7.3;. 

The  SAS  system  can  be  used  to  construct  confidence  bands  about  the  difference  of  two  models. 

Since  it  is  unlikely  that  the  values  of  the  covariates  will  provide  a  uniform  coverage  of  the  covariate 

region,  one  must  add  to  the  data  set  additional  observations  for  each  treatment  corresponding  to  the  desired 

values  of  the  covariate  (y  is  assigned  a  missing  value  for  these  observations).  The  following  SAS  system 

code  generates  one  set  of  values  to  be  added  to  the  data  set;  fits  three  models  and  constructs  the 

confidence  band  about  the  difference  between  the  models  of  treatments  1  and  2. 

DATA  RAW;  INPUT  TRT  y  x; 

CARDS 

(THE  DATA) 

DATA  GENERATE;  Y  =  ;; 
DO  TRT  =  1  TO  3; 

DO  X  =  1  TO  10; 

OUTPUT; 

END; 

END; 

DATA  ALL;  SET  RAW  GENERATE; 

PROC  GLM;  CLASS  TRT; 

MODEL  Y  =  TRT  X*TRT/SOLUTION; 

OUTPUT  OUT  =  VALUE  P  =  Y  STDP  =  STD; 

DATA  ONE;  SET  VALUE;  IF  TRT  =  1; 

PYl  =  PY;  STDl  =  STD; 

PROC  SORT;  BY  X; 

DATA  TWO;  SET  VALUE;  IF  TRT  =  2; 

PY2  =  PY;  STD2  =  STD; 

PROC  SORT;  BY  X; 

DATA  ONETWO;  MERGE  ONE  TWO;  BY  X; 
Q  =  SQRT(2*INVF(2,27)); 
DIF  =  PYl  -  PY2; 

STEDIF  =  SQRT(STD1**2  +  STD2**2); 
LOW  =  DIFF  -  Q*STEDIF; 
HIGH  =  DIFF  +  Q*STEDIF; 
PROC  PLOT;  PLOT 

DIFF*X  =  '*'  LOW*X  =  '-'  HIGH*X  =  '-' 
/OVERLAY; 

The  value  of  Q  in  the  SAS  system  code  corresponds  to  the  Scheffe  percentage  point.  The  analysis  and 
the  SAS  System  computations  described  in  this  section  are  demonstrated  in  detail  by  examples  in  the  next 
sections. 


2.8      Example:  Sugar  Beets  and  Nitrogen-Equal  Slopes 

In  a  study  of  growth  regulators  for  sugar  beets,  it  was  determined  that  there  was  substantial  plot 
to  plot  variation  in  the  level  of  available  nitrogen.  The  amount  of  nitrogen  in  the  soil  can  affect  the  yield 
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of  the  beets  in  addition  to  an  effect  due  to  tlie  growth  regulators.  After  planting  the  sugar  beets,  the 
available  nitrogen  was  measured  from  soil  samples  obtained  from  each  plot.  Since  there  is  a  lot  of  plot 
to  plot  variation  in  the  amount  of  available  nitrogen,  the  experimenter  wanted  to  be  able  to  compare  the 
effect  of  the  growth  regulators  at  specified  levels  of  nitrogen.  Thus,  the  available  nitrogen  in  the  soil  is 
used  as  the  covariate  and  the  dependent  variable  is  the  yield  of  the  sugarbeet  roots  per  plot  in  pounds. 
The  experimental  design  is  a  one-way  treatment  structure  with  three  growth  regulators  (TRT)  and  10  plots 
per  treatment  in  a  completely  randomized  design  structure.  The  data  are  in  Table  2.8.1. 


Table  2.8.1.  Plot  Yield  and  Residual  Nitrogen  for  the  Sugar  Beet  Example 


TRT  1  2  3 


YIELD 

NITROGEN 

YIELD 

NITROGEN 

YIELD 

NITROGEN 

210 

100 

155 

65 

155 

55 

150 

50 

150 

55 

160 

55 

200 

105 

170 

65 

165 

65 

180 

75 

175 

75 

185 

70 

190 

80 

185 

75 

185 

80 

220 

110 

195 

90 

200 

85 

170 

60 

205 

84 

205 

90 

170 

70 

210 

100 

215 

105 

190 

90 

215 

95 

220 

105 

220 

100 

220 

100 

220 

100 

The  SAS  system  statements  required  to  carry  out  the  analysis  are  given  in  Table  2.8.2  along  with 
an  annotation  of  each  statement's  function.  The  analysis  of  variance  tables  obtained  and  the  estimates  of 
the  parameters  are  given  in  Tables  2.8.3  to  2.8.5. 


Table  2.8.2.  SAS  System  Statements  for  Sugar  Beet  -  Nitrogen  Example 

i)  Read  in  the  data 

DATA  GOV;  INPUT  TRT  YIELD  NIT; 
CARDS; 

ii)  Fit  the  model  to  test  H    in  (2.6) 

PROG  GLM;  GLASSES  TRT; 

MODEL  YIELD=TRT  NIT*TRT  /  NOINT  SOLUTION; 

iii)  Fit  model  to  test  H    in  (2.11) 

PROG  GLM;  GLASSES  TRT; 
MODEL  YIELD=TRT  NIT  NIT*TRT; 

iv)  Fit  parallel  lines  model  to  compare  adjusted  means 

PROG  GLM;  GLASSES  TRT; 

MODEL  YIELD=TRT  NIT  /  SOLUTION; 

LSMEANS  TRT  /  STDERR  PDIFF; 
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Table  2.8J      Results  from  Fitting  Model  [2.1]  to  Data  for  Example  1  for  Parameter 
Estimation  and  to  Test  Ho,  of  (2.6).  From  Part  (ii)  of  Table  2.8.2. 

DEPENDENT  VARIABLE:  YIELD 


SOURCE 

MODEL 
ERROR 


DP 
6 

24 


SUM  OF  SQUARES 

1094097.45 
1101.55 


MEAN  SQUARE 

182349.56 
45.64 


SOURCE 


DP   TYPEIII  SS 


F VALUE 


PR>F 


TRT 
NIT*TRT 


1079250.00 
14847.45 


7830.94 
107.73 


.0001 
.0001 


PARAMETER 
TRT 


ESTIMATE 

98.943 
66.320 
88.780 


STD  ERROR  OF  EST 

9.61 
11.64 
9.74 


NIT*TRT 


1.084 
1.512 
1.262 


0.11 
0.14 
0.12 


Table  2.8.4.     Results  of  Fitting  Model  (2.1)  to  Data  for  Sugar  Beet  -  Nitrogen  Example  for 
Parallel  Line  Analysis 

(MODEL  and  ERROR  sums  of  squares  and  mean  squares  are  the  same  as  those  given  in  Table 
2.8.3.) 

SOURCE         d.f.  TYPE  IV  SS       F  VALUE  PR>F 

TRT  2  217.278  2.36  0.1155 

NIT  1  14726.736        320.57  0.0001 

NIT*TRT        2  257.555  2.80  0.0805 
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Table  2.8^      Results  from  Fitting  Model  (2.13)  to  Data  for  Sugar  Beet  -  Nitrogen  Example  for 
Parallel  Line  Analysis 


DEPENDENT  VARIABLE:  YIELD 
SOURCE  DP 


SUM  OF  SQUARES 


MODEL 
ERROR 
CORRECTED 
TOTAL 


SOURCE 
TRT 


3 

26 
29 

DP 

2 
1 


14636.56 
1360.10 
15996.66 


PARAMETE 
R 

INTERCEPT 
TRT  1 
TRT  2 
TRT  3 
NIT 


TRT 


ESTIMATE 


89.56  B 
-4.76  B 
-2.37  B 

0.00  B 

1.25 


OBSERVED 
MEAN 


T  FOR  HO: 
PARAMETER=0 


13.80 
-1.47 
-0.73 

16.70 

ADJUSTED 
MEAN 


TYPEIV  SS 

112.60 
14589.89 


F VALUE 

1.08 
278.90 


MEAN 
SQUARE 

4878.85 
52.31 


PR>F 

0.3556 
0.0001 

STD  ERROR 
OF  ESTIMATE 


6.49 
3.24 
3.23 

0.07 

PROB> I T I 
HO:LSMEAN=0 


pr>|t| 

0.0001 
0.1543 
0.4696 

0.0001 
STDERR(ADJ  MN) 


1  190.0  187.28                         2.293  0.0001 

2  188.0  189.67                        2.289  0.0001 

3  191.0  192.04                        2.288  0.0001 

t-tests  for  differences  between  adjusted  mean 
PROB  >  I  T  I  HO:  LSMEAN(I)  =  LSMEAN(J) 

I/J  1               2  3 

1  .  0.4693  0.1543 

2  0.4693  .  0.4696 

3  0.1543  0.4696 


Model  (2.1)  was  fit  to  the  data  with  the  results  appearing  in  Table  2.8.3.  The  sums  of  squares  for 
NIT*TRT  tests  Hqi  of  (2.6)  and  is  highly  significant.  Thus  the  experimenter  would  conclude  that  the 
expected  mean  yields  do  depend  on  NIT  (available  nitrogen).  Table  2.8.3  also  displays  die  estimates  of 

the  model  parameters,     's  correspond  to  TRT  and  l^^'s  correspond  to  NIT*TRT. 

Next,  the  equality  of  slopes  hypothesis  H02  of  (2. 1 1)  is  tested  to  determine  if  a  parallel  lines  model 
will  adequately  describe  the  data.  The  sum  of  squares  due  to  NIT*TRT  in  Table  2.8.4  is  SSH02.  The 
significance  level  of  the  test  statistic  is  .0805.  If  the  decision  is  made  at  the  .05  level,  then  one  would 


37 


conclude  that  a  parallel  lines  model  is  adequate  to  describe  the  data.  However,  if  the  decision  is  made 
at  the  .10  level,  then  one  would  conclude  that  the  parallel  lines  model  is  not  adequate.  For  the  purpose 
of  discussion,  the  .05  level  is  used  and  the  parallel  lines  model  is  determined  to  be  adequate  to  describe 
the  relationship  between  yield  and  level  of  available  nitrogen  (one  should  plot  the  residuals  before 
continuing). 

The  parallel  lines  model  of  (2.13)  was  fit  to  the  data  and  the  results  are  given  in  Table  2.8.5. 

The  sum  of  squares  due  to  TRT  is  testing  H03  of  (2.17)  given  that  H02  is  true.  The  significance 
level  is  .3556,  indicating  that  the  distances  between  the  regression  lines  are  not  significantly  different  from 
zero.  The  estimated  slope  of  the  parallel  regression  lines  is  1.252,  which  is  significantly  different  from 
zero. 

Table  2.8.5  also  contains  the  observed  TRT  means  for  yield  and  the  adjusted  means  of  the  TRT 
yields.  The  adjusted  means,  provided  by  the  LSMEANS  statement,  are  the  predicted  values  obtained  from 
the  treatment  regression  lines  evaluated  at  the  mean  NIT  (nitrogen)  value  for  the  experiment.  The  mean 
NIT  value  for  the  experiment  is  81.83.  For  example,  the  adjusted  mean  for  TRT  1  is 

Myjnit  -  81.833  =  (89.56  -  4.76)  +  81.83(1.25)  =  187.65. 
The  quantity  (89.56  -  4.76)  =  84.80  is  the  intercept  (ftj)  for  treatment  1  from  the  parallel  lines  model. 

The  significance  levels  of  the  t-tests  for  comparing  each  pair  of  adjusted  treatment  means  is 
presented  in  Table  2.8.5.  The  significance  levels  are  obtained  by  including  the  PDIFF  option  in  the 
LSMEANS  statement.  The  t-statistics  to  compare  adjusted  means  for  treatments  i  and  i  are  computed  as 

^  ^  Difference  between  two  adjusted  means 

Standard  error  of  the  difference  of  two  adjusted  means 

^^Y.lNIT  «  81.833  ~  ^^Y,|NrT  -  81.833 
^(l^.lNIT  -  81.833  -  AtJNIT  -  81.83333 

where  the  standard  error  is  computed  via  (2.23). 

The  probability  values  can  be  used  to  construct  an  LSD  procedure  by  concluding  the  distance 
between  two  lines  is  significantly  different  from  zero  if  the  printed  significance  level  is  less  than,  say,  .05. 
A  Bonferroni  multiple  comparison  procedure  can  be  used  by  concluding  the  distance  between  two  lines 
is  significantly  different  from  zero  if  the  printed  significance  level  is  less  than  .05/3  =  .01667.  Of  course, 
the  above  decisions  could  be  made  at  some  significance  level  other  than  .05. 

For  this  data  set,  none  of  the  lines  are  significantly  different  from  each  other,  therefore,  one  line  could  be 
used  to  describe  all  of  the  data.  Thus,  conclude  that  there  are  no  differences  between  the  growth  regulator 
treatments.  A  graph  of  the  data  and  estimated  regression  lines  are  in  Figure  2.4. 
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2.9      Example-Exercise  Programs  and  Initial  Resting  Heart  Rate 

An  exercise  physiologist  structured  three  types  of  exercise  programs  (EPRO)  and  conducted  an 
experiment  to  evaluate  and  compare  the  effectiveness  of  each  program.  The  experiment  consisted  of 
subjecting  a  person  to  a  given  exercise  program  for  eight  weeks.  At  the  end  of  the  training  program,  each 
person  ran  for  six  minutes  after  which  their  heart  rate  was  measured.  An  exercise  program  is  determined 
to  be  more  effective  if  individual's  on  that  program  have  lower  heart  rates  than  individuals  on  another 
exercise  program.  Since  people  entered  the  experiment  at  differing  degrees  of  fitness,  the  resting  heart 
rate  before  beginning  training  was  recorded,  and  was  used  as  a  covariate.  The  object  of  the  study  is  to 
be  able  to  compare  exercise  programs  at  a  common  initial  resting  heart  rate.  To  carry  out  the  experiment, 
24  males  between  28  and  35  years  of  age  were  selected  and  then  8  males  were  randomly  assigned  to  each 
of  the  three  EPRO  treatments.  The  exercise  program  (EPRO),  heart  rate  (HR)  after  the  6  minutes  run  at 
the  completion  of  eight  weeks  of  training  and  the  initial  resting  heart  rate  (IHR)  are  given  in  Table  2.9.1. 


Table  2.9.1  Data  for  Exercise  Program  Example 

EXERCISE  PROGRAM 
1  2  3 


HR 

IHR 

HR 

IHR 

HR 

IHR 

118 

56 

148 

60 

153 

56 

138 

59 

159 

62 

150 

58 

142 

62 

162 

65 

158 

61 

147 

68 

157 

66 

152 

64 

160 

71 

169 

73 

160 

72 

166 

76 

164 

75 

154 

75 

165 

83 

179 

84 

155 

82 

171 

87 

177 

88 

164 

86 

The  SAS  system  statements  necessary  to  analyze  this  data  set  are  given  in  Table  2.9.2.  Model 
2.1  was  first  fit  to  the  data  and  the  results  appear  in  Table  2.9.3.  The  sum  of  squares  for  the  IHR*EPRO 
tests  Hoi  of  (2-6)  and  is  highly  significant.  Thus,  conclude  that  the  heart  rate  depends  on  the  initial  resting 
hear  rate.  Table  2.9.3  also  contains  the  estimates  of  the  parameters, 

fiCj's  corresponding  to  EPRO  and  ^.'s  corresponding  to  IHR*EPRO. 
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Table  2.9.2  SAS  System  Statements  for  Exercise  Program  Example 

i)  Read  in  the  data  cards. 

DATA  COVI;  INPUT  EPRO  HR  IHR;  CARDS; 

ii)  Fit  model  to  test  Hqi  in  (2.6)  and  to  provide  estimate  of  several  linear  combinations  of  the 
parameters 

PROC  GLM:  CLASSES  EPRO; 

MODEL  HR=EPRO  IHR*EPRO/NOINT  SOLUTION 

LSMEANS  EPRO/STDERR  PDIFF; 

ESTIMATE  'B1-B2,  IHR*EPRO  1  -1  0; 

ESTIMATE  'B1-B3'  IHR*EPRO  1  0  -1; 

ESTIMATE  'B2-B3'  IHR*EPRO  0  1-1: 

ESTIMATE  '1  AT  60  -  2  AT  80;  EPRO  1-10  IHR*EPRO  60  -80  0: 
ESTIMATE  'A1-A2'  EPRO  1  -1  0; 
ESTIMATE  'A2-A3'  EPRO  0  1-1; 
ESTIMATE  'A1-A3'  EPRO  10-1; 

ESTIMATE  '1-2  AT  55'  EPRO  1-10  IHR*EPRO  55  -55  0; 
ESTIMATE  '1-2  AT  70'  EPRO  1  -1  0  IHR*EPRO  70  -70  0; 
ESTIMATE  '1-2  AT  85'  EPRO  1-10  IHR*EPRO  85  -85  0; 
ESTIMATE  '1-3  AT  55'  EPRO  1  0  -1  IHR*EPRO  55  0  -55; 
ESTIMATE  '1-3  AT  70'  EPRO  10-1  IHR*EPRO  70  0  -70; 
ESTIMATE  '1-3  AT  85'  EPRO  1  0  -1  IHR*EPRO  85  0  -85; 
ESTIMATE  '2-3  AT  55'  EPRO  0  1  -1  IHR*EPRO  0  55  -55; 
ESTIMATE  '2-3  AT  70'  EPRO  0  1-1  IHR*EPRO  0  70  -70; 
ESTIMATE  '2-3  AT  85'  EPRO  0  1  -1  IHR*EPRO  0  85  -85; 

iii)  Fit  model  to  test      in  (2.1 1) 

PROC  GLM;  CLASSES  EPRO; 

MODEL  HR  =  EPRO  IHR  IHR*EPRO  /  SOLUTION  NOINT; 


Table  2.93      Results  From  Fitting  Model  (2.1)  to  the  Data  in  the  Exercise  Program  Example 
for  Parameter  Estimation  and  to  Test  Hoi  of  (2.6)  [Part  (ii)  of  Table  2.9.2] 

DEPENDENT  VARIABLE:  HR 

SUM  OF 


SOURCE 

DF 

SQUARES 

MEAN  SQUARE 

MODEL 

6 

594974.35 

99162.38 

28.20 

ERROR 

18 

507.64 

UlNL-UKKilL.  1 E/U 

ERROR 

SOURCE 

DF 

TYPE  IV  SS 

F  VALUE 

PR>F 

EPRO 

3 

5108.68 

60.38 

0.0001 

IHR*EPRO 

3 

2650.60 

31.33 

0.0001 

T  FOR  HO: 

oTD  ERROR 

PARAMETER 

ESTIMATE 

PARAMETER=0 

PR>|T| 

OF  ESTIMATE 

EPRO 

1 

46.53 

3.66 

0.0018 

12.71 

2 

97.19 

6.88 

0.0001 

14.12 

3 

137.48 

10.97 

0.0001 

12.52 

IHR*EPRO 

1 

1.48 

8.29 

0.0001 

0.17 

2 

0.93 

4.80 

0.0001 

0.19 

3 

0.26 

1.47 

0.1576 

0.17 

Next  the  experimenter  needs  to  determine  if  a  parallel  lines  model  will  adequately  describe  the 
data,  i.e.,  test  H02  of  (2.11).  The  sum  of  squares  due  to  IHR*EPRO  in  Table  2.9.4  is  SSH02  (from  part 
iii  of  Table  2.9.2).  The  significance  level  of  the  test  is  .0(X)6,  indicating  that  the  parallel  lines  model  is 
not  reasonable  to  describe  this  data.  Thus,  it  is  necessary  to  compare  the  EPRO  at  various  values  of  IHR 
by  using  the  unequal  slope  model.  The  experimenter  chose  to  compare  the  three  exercise  programs  at  IHR 
=  55,  70  and  85.  In  Table  2.9.2,  the  last  9  ESTIMATE  statements  are  those  required  to  make  the  three 
pairwise  comparisons  between  the  EPR's  at  each  given  value  of  IHR.  For  example,  '1-2  AT  55'  asks 
GLM  to  compare  EPRO  1  with  EPRO  2  at  IHR  =  55,  etc.  The  results  of  these  nine  tests  are  in  Table 
2.9.5.  Table  2.9.5  also  contains  the  results  of  testing  the  following  special  hypotheses: 

i)        Pi  =  P2'  Pi  =  Ps  and  p2  =p3 
iii)       a,  =  02',  ttj  =  (X3  and  0^  =  0^ 
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Investigating  the  relationship  between  the  slopes  is  reasonable  as  one  might  wish  to  simplify  the  model 
by  grouping  treatments  with  "equal"  slopes  into  groups  where  a  parallel  lines  model  is  fit  to  the  data 
within  a  group  and  nonparallel  lines  between  groups.  When  such  a  simplification  can  occur,  the 
comparison  process  is  easier  as  one  can  compare  lines  within  a  group  with  the  LSMEANS. 
Comparing  two  treatments  at  different  values  of  the  covariate  (ii  from  above)  is  not  usually  done,  but  there 
are  some  circumstances  where  such  comparisons  are  warranted  Gike  comparing  EPRO  1  at  60  to  EPRO  2 
at  80,  not  meaningful  here).  But,  for  example,  suppose  it  costs  die  same  for  Xi  units  of  the  covariate  on 
process  1  as  for  X2  units  of  the  covariate  on  process  2,  then  it  is  reasonable  to  compare  the  models  at 
different  values  of  X,  but  on  an  equal  cost  basis.  Comparison  of  the  intercepts  (part  iii  from  above) 
compares  the  distances  between  the  regression  lines  at  IHR  =  0.  In  this  example,  it  is  stupid  to  compare 
the  intercepts.  We  must  be  careful  not  be  make  such  stupid  comparisons. 

Table  2.9.4      Results  of  Fitting  Model  (2.1)  to  Exercise  Program  Example  to  Test  H02  of  (2.6) 

(Model  and  ERROR  sum  of  squares  and  the  mean  squares  are  the  same  as  Table  2.8. 

SOURCE       DF  TYPE  IV  SS  F  VALUE  PR>F 

EPRO  3  3108.68       60.38  0.0001 

IHR  1  1990.74        70.59  0.0001 

IHR*EPRO     1  658.98        11.68  0.0006 


The  adjusted  means  in  Table  2.9.5  estimate  the  EPRO  means  at  IHR  =  70.375  the  overall  mean 
IHR  for  the  experiment.  While  it  might  be  of  interest  to  compare  the  treatment  means  at  70.375,  they 
should  also  be  compared  to  other  values  of  the  covariate  as  well.  Figure  2.5  summarizes  the  analysis  with 
the  estimated  regression  lines  and  the  comparisons  of  means  at  IHR  =  55,  70  and  85.  If  two  means  at  the 
same  value  of  IHR  are  followed  by  the  same  letter,  the  means  are  not  significantly  different  while  unlike 
letters  indicate  a  difference.  A  .05  Bonferroni  approach  was  used  at  each  level  of  IHR  where  a  difference 
is  declared  significant  within  a  IHR  value  if  the  t-test  has  a  significance  level  less  that  or  equal  to  .05/3 
=  .01667. 

Another  way  to  compare  the  regression  lines  is  to  construct  confidence  bands  about  the  difference 
of  each  pair  of  models  and  determine  the  values  of  the  independent  variable  where  the  bands  exclude  zero. 
To  provide  control  over  the  error  rates,  Scheffe  percentage  points  are  used  in  constructing  the  bands.  The 
difference  between  two  models  (ith  and  i'th)  is  (a^  -  ct^)  +  (p,  -  P2)X  which  depends  on  two  parameters 
(tti  -  0(2)  and  (pi  -  P2).  Thus  the  percentage  point  used  is 

(2 
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Plot  of  HR*IHR.       Symbol  used  is  'o'. 
Plot  of  PHR*IHR.     Symbol  is  value  of  EPRO. 
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Figure  2.5  Graphs  of  the  estimated  models  for  the  exercise  program  example.  Bonferroni 
comparisons  between  the  exercise  programs  are  shown  for  IHR  =  55,  70,  85,  where 
models  with  the  same  letter  are  not  significantly  different  at  (x=.05.  Comparisons  are 
obtained  from  Table  2.9.5. 


43 


Table  2.9.5  Estimates  of  Selected  Linear  Combinations  of  the  Parameters  for  Exercise  Program  Example 


T  v:f\o  vjn- 
1  rL>K  nu: 

o  1 U  tiKKUK 

PARA  MPTFR 

Co  1  U.Vl/\  1  c 

PAR  AMFTFR— n 

PR>  1  T  1 

B1-B2 

0.54 

2.06 

0.0537 

0.26 

B1-B3 

1.22 

4.83 

0.0001 

0.25 

B2-B3 

0.67 

2.54 

0.0203 

0.26 

1  AT  60  -  2  AT  80 

-36.57 

-10.11 

0.0001 

3.61 

A1-A2 

-50.65 

-2.66 

0.0158 

19.01 

A2-A3 

-40.29 

-2.13 

0.0469 

18.88 

A1-A3 

-90.94 

-5.09 

0.0001 

17.85 

1-2  AT  55 

-20.55 

-4.11 

0.0007 

5.00 

1-2  AT  70 

-12.34 

-4.62 

0.0002 

2.67 

1-2  AT  85 

-4.13 

-0.91 

0.3769 

4.56 

1-3  AT  55 

-23.76 

-5.19 

0.0001 

4.58 

1-3  AT  70 

-5.44 

-2.05 

0.0555 

2.65 

1-3  AT  85 

12.87 

2.75 

0.0132 

4.68 

o  q  AT 

-kj.oj 

2-3  AT  70 

6.90 

2.58 

0.0190 

2.67 

2-3  AT  85 

17.01 

3.64 

0.0019 

4.67 

HR 

STD  ERR         PROB> | T | 

LSMEAN 

EPRO 

LSMEAN 

LSMEAN  HO:LSMEAN=0 

NUMBER 

1 

151.06 

1.87  0.0001 

1 

2 

163.20 

1.89  0.0001 

2 

3  156.04  1.88  0.0001  3 


PROB  >  I  T I  HO:  LSMEAN(I)=LSMEAN(J) 
yj            1  2  3 

1  0.0002  0.0775 

2  0.0002  0.0154 

3  0.0775  0.0154 


where  x>  is  the  degrees  of  freedom  associated  with  the  MSERROR.  Since  the  estimates  of  the  parameters 
for  the  models  are  independent  between  the  models,  the  variance  of  the  difference  of  two  models  is  the 
sum  of  the  variances  of  the  two  predicted  values.  Most  computer  codes  provide  the  predicted  values  and 
the  standard  errors  of  the  predicted  values  which  can  be  combined  to  construct  the  confidence  band  about 
the  difference  of  two  models.  The  end  points  of  the  confidence  band  at  a  given  value  of  the  independent 
variable  is 

(d,-H|i^X)  -  (6^^^^,X)  ±  (2  F„^^p>ar(d,+3,X)  -h  var((3^-H(5,X))^ 

The  annotated  SAS  system  code  used  to  generate  the  data  sets  at  values  of  IHR  from  55  to  85 
(and  generate  the  necessary  missing  values  for  HR)  and  construct  the  confidence  bands  is  displayed  in 
Table  2.9.6.  The  three  graphs  for  comparing  all  pairs  of  models  are  in  Figures  2.6  to  2.8. 
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Table  2.9.6      SAS  System  code  for  constructing  confidence  bands  about  the  difference  of  pairs  of 
treatment  (EPRO)  models. 

a)  Generate  a  data  set  with  IHR  values  from  55  to  80  for  each  EPRO. 
DATA  PLOT;  HR=.; 

DO  EPRO=  1  TO  3; 

DO  IHR=55  TO  85; 

OUTPUT; 

END; 

END; 

b)  Merge  the  generated  data  with  the  actual  data  set. 
DATA  ALL;  SET  HEART  PLOT; 

c)  Use  PROC  GLM  to  estimate  the  parameters  of  the  models  and  compute  the  predicted  values  and 
their  standard  errors.  Write  the  results  to  a  data  set  --  XPLOT. 

PROC  GLM;  CLASSES  EPRO; 

MODEL  HR=EPRO  IHR*EPRO/NOINT  SOLUTION; 

OUTPUT  OUT=XPLOT  P=PHR  STDP=STDR; 

d)  For  data  sets  for  each  treatments  data. 

DATA  ONE;  SET  XPLOT;  IF  EPR0=1  AND  HR=.;P1=PHR;S1=STDR; 
PROC  SORT;  BY  IHR; 

DATA  TWO;  SET  XPLOT;  IF  EPR0=2  AND  HR=.;P2=PHR;S2=STDR; 
PROC  SORT;  BY  IHR; 

DATA  THR;  SET  XPLOT;  IF  EPR0=3  AND  HR=.;P3=PHR;S3=STDR; 
PROC  SORT;  BY  IHR; 

e)  Merge  pairs  of  data  sets  by  the  IHR,  compute  the  difference  between  the  models  at  each  IHR, 
compute  the  standard  error  of  each  difference,  compute  the  upper  and  lower  limits,  and  plot  the 
results.  The  value  2.665  =(2  *  3.55)^,  the  Scheffe  percentage  point  for  each  pair  of  models.  The 
SAS  system  code  is  included  for  comparing  models  1  and  2. 

DATA  ONE_TWO;  MERGE  ONE  TWO;  BY  IHR; 
DIFF=P1  -P2;STE=SQRT(S  1  *S  1+S2*S2);  Z=0; 
LOW=DIFF-2.665*STE; 
HIGH=DIFF+2.66*STE; 

PROC  PLOT;  PLOT  DIFF*IHR='*'  LOW*IHR='+'  HIGH*IHR='+'  Z*IHR='-'  /OVERLAY; 
TITLE3  'TREATMENT  1  MINUS  TREATMENT  2'; 
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Plot  of  DIFF*IHR.     Symbol  used  is 
Plot  of  LOW*IHR.       Symbol  used  is  '+'. 
Plot  of  HIGH*IHR.     Symbol  used  is   ' +' . 
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Figure  2.6       Confidence  band  about  the  difference  of  the  models  for  EPRO  one  and  two. 
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Plot  of  DIFF*IHR.     Symbol  used  is 
Plot  of  LOW^IHR.       Symbol  used  is 
Plot  of  HIGH*IHR.     Symbol  used  is   ' +' . 
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Figure  2.7       Confidence  band  about  the  difference  of  the  models  for  EPROs  one  and  three. 
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Plot  of  DIFF^IHR. 
Plot  of  LOW^IHR. 
Plot  of  HIGH*IHR. 
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Figure  2.8 


Confidence  band  about  the  difference  of  the  models  for  EPROs  two  and  three. 
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Summary 

One  must  not  blindly  use  computer  code  to  do  analysis  of  covariance.  Several  issues  which  must 
be  addressed  are  (1)  make  sure  the  model  adequately  describes  the  data  Ginear,  etc.),  (2)  check  the  equality 
of  variances,  (3)  check  the  equality  of  slopes,  (4)  use  the  right  model  to  make  comparisons  between  the 
treatments,  and  (5)  make  sure  you  know  what  information  is  being  obtained  from  your  computer  code. 
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Abstract 

Solutions  to  two-phase  regression  problems  are  generally  found  by  using  either  the  conventional 
method,  which  is  a  series  of  linear  regressions,  or  alternately  by  using  a  nonlinear  regression  approach. 
The  conventional  method  even  in  the  simplest  cases  is  time  consuming  prompting  the  researcher  to  turn 
to  a  quick  solution  via  nonlinear  regression.  The  use  of  one-shot  nonlinear  regression  to  solve  two-phase 
regression  problems  often  leads  to  incorrect  solutions  that  look  reasonable.  To  ensure  that  nonlinear 
regression  reaches  the  correct  solution  bounds  must  be  placed  on  one  of  the  predictors. 

Introduction 

Two-phase  straight-line  regression  is  a  statistical  tool  that  has  wide  applications  in  many  fields 
of  research,  and  particularly  in  the  areas  of  plant  and  animal  physiology.  A  recent  study  (Nickerson, 
Facey,  Grossman  1989)  indicates  that  the  two-phase  procedure  is  only  seldom  used  to  its  full  potential  on 
suitable  data. 

In  this  paper  I  will  review  the  Nickerson  et  al.  fmdings.  Further,  1  will  summarize  the 
conventional  least-squares  approach  to  two-phase  regression.  Finally,  I  will  discuss  a  modified  nonlinear 
approach  and  give  an  example  using  the  SAS  procedure  NLIN. 

How  Well  is  Two-Phase  Regression  Being  Used? 

In  order  to  obtain  an  estimate  of  how  appropriately  two-phase  regression  is  being  used  over  a  wide 
range  of  data  Nickerson,  Facey  and  Grossman  (1989)  studied  all  the  papers  published  in  Physiological 
Zoology  from  1983  to  1987.  They  found  that  of  the  31  papers  that  explicitly  dealt  with  data  suited  to 
two-phase  regression  only  one  paper  provided  a  statistically  valid  description  of  its  data  set. 

The  remaining  30  papers  were  categorized  as  follows:  eight  had  no  precise  definition;  nine  fitted 
a  regression  line  through  only  one  section  of  the  data;  five  did  a  separate  regression  of  each  section  of  the 
data  but  did  not  constrain  the  lines  to  meet  between  the  two  sets;  seven  did  a  visual  estimate  of  the  join 
point,  (the  join  point  is  defined  as  the  location  of  the  intersection  of  the  two  straight  lines  which  constitute 
the  two-phase  regression);  and  one  researcher  did  a  series  of  regressions  based  on  preassigned  x-values 
for  join  points  and  chose  the  combination  which  produced  the  smallest  sum  of  squares.  While  some  of 
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the  above  methods  seem  to  produce  fairly  reasonable  estimates  of  the  join  points  none  of  these  produces 
a  truly  least-squares  best  fit  estimate. 

Conventional  Method  of  Fitting  a  Two-Phase  Straight-line  Regression 

Specification  of  the  Model.  This  section  uses  equations  and  examples  from  Hudson  (1966)  and 
Nickerson,  Facey  and  Grossman  (1989).  Please  see  these  authors  for  a  fuller  description  of  the  topic.  For 
a  discussion  of  the  power  of  two-phase  regression  see  Hinkley  (1971). 

Given  a  set  of  data  as  in  Figure  1  and  ignoring  the  distinction  between  sample  and  population 
parameters,  the  equations  which  state  the  conventional  two-phase  regression  are 
yi  =  INTl  +  SLPl  X  Xi  -H  Ci,        if  Xj  <  XJOESF,  and 

(1) 

yi  =  INT2  +  SLP2  x  Xj  +  Cj,        if  x^  >  XJOIN. 
In  the  best  fit  solution  S  Q\  is  a  minimum.  Since  the  two  regression  lines  meet  at  the  join  point  whose 
abscissa  is  XJOIN,  it  follows  that 

(2)  INTl  +  SLPl  X  XJOIN  =  INT2  +  SLP2  x  XJOIN, 

and  the  parameter  of  XJOIN  can  be  calculated  directly  from  equation  (2), 

(3)  XJOIN  =  (E^^Tl  -  INT2)  /  (SLP2  -  SLPl). 


Fig.  1 
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Procedure  for  Determining  the  Best  Two-Phase  Model.  First  we  order  the  x's  from  lowest  to 
highest  so  that  Xj  <  X2  <  ...  <  x^.  We  next  split  the  data  into  two  complementary  sets  corresponding  to 
the  two  phases  in  which  one  set  consists  of  the  data  points  whose  abscissas  are  Xi  to  Xj  while  the  other 
set  has  abscissas  Xj+i  to  x^.  Simple  linear  regression  is  performed  on  each  of  the  two  sets.  Since  the 
minimum  number  of  observations  required  to  define  a  regression  line  is  two  then  2  <  j  <  n-2,  where  n  is 
the  number  of  observations.  The  maximum  number  of  values  for  j  is  n-3  and  is  reduced  by  one  for  each 
duplicate  x-value. 

The  coefficients  from  each  of  the  n-3  (fewer  if  there  are  x-value  duplicates)  pairs  of  regressions 
are  used  to  compute  the  corresponding  XJOIN  as  in  equation  (3).  How  we  proceed  from  here  depends 
on  whether  or  not  Xj  <  XJOINj  <  Xj+i.  If  XJOIN  for  any  j  lies  between  the  neighboring  abscissas  of  its 
disjoint  data  sets  then  we  have  an  ordinary  least-squares  (OLS)  solution.  Tlie  coefficients  for  such  a  j  are 
as  stated  in  equation  (1)  and  the  SSE  of  the  combined  equation  is  the  sum  of  the  SSEs  of  each  of  the  two 
phases.  If,  mirabile  dictu,  all  n-3  pairs  of  regressions  produce  OLS  solutions  then  the  best  solution 
belongs  to  the  pair  with  the  smallest  SSE. 

For  any  value  of  j  whose  XJOIN  does  not  lie  within  the  exclusive  range  of  Xj  to  Xj+i  the  solution 
is  not  valid.  Here  the  regression  lines  have  been  generated  with  one  or  more  data  points  placed  in  the 
wrong  phase.  This  invalid  pair  of  regression  lines  can  be  redefined  by  forcing  the  join  point  to  equal  Xj. 
This  is  called  a  constrained  least-squares  (CLS)  solution.  The  equations  for  transforming  an  invalid  OLS 
solution  into  a  CLS  solution  are  taken  from  Nickerson,  et  al.  In  the  following  definifions  n^  and  SS^  refer 
respectively  to  the  number  of  observations  and  sum  of  squares  in  the  first  phase  while  n^.  and  SS^'  refer 
similarly  to  the  second  phase. 

(4)  CLSINTl=INTl-(s/t)x(n^xSS^J-^xlJi^i[(Xi-Xj)xXi], 

(5)  CLSSLPl=SLPl+(s/t)x(njXSSJJ-^xSi^i(Xi-Xj), 

(6)  CLSINT2=INT2-(s/t)x(nj.xSSJ'J-^xri^^,[(Xi-Xj)xXi],  and 

(7)  CLSSLP2=SLP2+(s/t)x(n^.xSSJ'j-^xI",^,i(xrXj), 
where 

s  =  (DSrri  +  SLPl  X  Xj)  -  (INT2  +  SLP2  x  x^, 

and 

t=(njXSS^J-^xI^,.i(x,-x.)^+(nj.xSS^'J-^xI"i^,i(xpxpl 
For  the  Q.S  solution  the  sum  of  squares  for  error  =  S"=i  (yryf-   For  any  j -value  the  SSE  of  the 
constrained  solution  is  greater  than  the  SSE  of  the  corresponding  non-constrained  solution  be  it  valid 
(OLS)  or  invalid.  It  is,  therefore,  not  necessary  to  perform  CLS  adjustments  for  any  invalid  solutions 
whose  SSEs  are  greater  than  the  smallest  SSE  among  the  valid  solutions. 

Even  the  simplest  OLS  solution  requires  up  to  (n-3)x2  simple  linear  regressions  to  obtain 
preliminary  results  and  tiiereafter  it  requires  either  a  great  deal  of  hand  calculation  or  some  non-trivial 


53 


programming  to  accumulate  and  process  the  various  regression  parameters.  In  the  more  typical  problem 
involving  some  CLS  solutions  the  time  required  to  do  the  calculations  could  run  into  many  hours  while 
the  risk  of  computational  errors  grows  with  the  complexity  of  the  problem.  In  those  experiments  for 
which  regressions  are  required  for  each  of  many  experimental  units  the  calculation  time  can  be  very 
lengthy. 

A  Nonlinear  Approach 

Nonlinear  procedures  such  as  those  contained  in  the  SAS,  SPSS-x,  and  BMDP  packages  are  well 
suited  to  finding  quick  and  accurate  solutions  to  problems  for  which  the  SSE  is  a  smooth  function  of  each 
of  the  predictors.  In  two-phase  regression  the  SSE  is  not  a  smooth  function  of  XJOIN.  A  modified 
nonlinear  procedure  is  therefore  required.  In  this  section  I  will  employ  the  SAS  procedure  NLIN  to 
illustrate  this  methodology.  This  procedure  requires  two  adjustments  not  found  in  the  conventional 
method:  first,  all  the  data  are  submitted  to  the  program  as  one  set  rather  than  being  run  twice  as  two 
disjoint  sets,  and  second,  since  the  SSE  as  a  function  of  XJOIN  is  not  smooth  at  values  XJOIN=Xi 
(Hudson  1966,  p.  1105),  it  is  necessary  to  constrain  the  regression  so  that  XJOIN  does  not  equal  any  Xj. 

In  SAS's  NLIN  procedure  it  is  possible  to  process  both  sets  of  data  as  defined  in  equation  (1) 
simultaneously  within  one  regression  by  using  the  IF  statement  to  branch  to  either  data  set  and  its 
appropriate  MODEL  and  DERivative  statements  (see  SAS  1990,  pp.  1 162-3).  For  the  sake  of  generality, 
I  will  redefine  equation  (1)  so  that  only  one  model  statement  need  be  given.  The  single  model  statement 
can  be  used  on  a  package  lacking  branching  capability  provided  it  has  a  derivative-free  nonlinear 
regression  routine. 

The  new  parameters  are  defined  as  follows:  YJOIN  is  the  ordinate  of  the  join  point;  XJOIN  is 
the  abscissa  of  the  join  point;  AVESLOPE  is  the  average  of  SLPl  and  SLP2  from  equation  (1);  and 
DELSLOPE  is  SLP2  minus  AVESLOPE.  The  model  may  then  be  restated  as 
(8)      yi  =  YJOIN  +  AVESLOPE  x  (x^  -  XJOIN)  +  DELSLOPE  x  ABS(Xi  -  XJOIN). 
In  the  above  equation  ABS  is  a  function  that  returns  the  absolute  value  of  the  expression  in  parentheses. 

Example 

The  SAS  program  listed  below  does  regression  separately  for  two  experimental  units.  Cow  A  and 
Cow  B.  The  data,  which  are  borrowed  from  Hudson  (1966),  permit  accurate  checking  of  results. 


Cow  A  y 
Cow  B  y 


X 

1 

2 

3 

4 

5 

1 

2 

4 

4 

3 

1 

1 

2 

4 

7 

3 

1 
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For  each  data  set  the  program  searches  for  a  minimum  SSE  with  XJOEN  constrained  to  lie 
successively  within  each  of  the  three  data  intervals.  The  intervals  are  defined  by  the  BOUNDS  statement. 
One  must  be  careful  here  to  shrink  the  range  between  the  x-value  limits  because  SAS  treats  the  hard 
inequality  signs,  <  and  >,  as  if  they  were  soft  inequality  signs,  <  and  >.  The  MODEL  statement  is  defined 
by  equation  (8). 

Pure  mathematicians  will  be  shocked  to  see  that  DELSLOPE  has  a  partial  derivative,  ABS(HOUR- 
XJOIN).  The  program  is  not  affected  by  the  lack  of  smoothness  in  the  derivative  because  the  point  at 
which  the  slope  changes  abruptly,  HOUR=XJOIN,  has  been  precluded  from  any  computations  by  the 
bounds  put  on  XJOIN. 

Table  1  contains  the  parameters  (regrouped  by  cow  ID)  that  were  written  by  the  PUT  statement 
at  the  end  of  each  of  the  three  intervals.  Comparing  the  results  of  Cow  A  with  the  curve  in  Figure  2a, 
one  notes  that  the  NLIN  procedure  identifies  the  minimum  in  the  second  interval  and  the  local  minimum 
in  the  third.  Both  are  OLS  solutions  and  are  accurate  to  about  eight  digits.  For  Cow  B  the  second  and 
third  interval  solutions  place  XJOIN  at  the  predefined  boundaries.  These  are  CLS  solutions.  The  size  of 
the  error  in  the  computed  value  of  each  of  the  parameters  is  a  function  of  the  size  of  the  buffer  around 
the  x-values  at  the  boundaries. 


Table  1 


Interval 

Cow 

XJOIN 

YJOIN 

SLOPE  1 

SLOPE  2 

SS  Error 

1 

A 

2.999990 

4.315794 

1.78949 

0.92105 

1.39475111 

2 

A 

3.611111 

4.750000 

1.50000 

-1.50000 

0.33333333 

3 

A 

4.193548 

4.612903 

1.10000 

-2.00000 

0.70000000 

1 

B 

2.999990 

5.368425 

2.42107 

-1.15789 

11.57898127 

2 

B 

3.999990 

6.315802 

1.92107 

-2.78946 

1.39477651 

3 

B 

4.000010 

6.315801 

1.92105 

-2.78950 

1.39473684 

Hudson 

4.000000 

6.315789 

1.92105 

-2.78947 

1.39473684 

Listing  of  SAS  Nonlinear  Program 

DATA  raw; 

INFILE  'a:twophase.dat'; 

INPUT  cow  $  hour  hormone  KLICK; 

PROC  SORT;  BY  cow; 
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/*     ****  Solve  for  first  interval    ****  */ 
PROC  NLIN  DATA=raw;  BY  cow; 

FARMS  YJOIN=  5  AVESLOPE=  0.1  DELSL0PE=-2      XJOIN=  2.5  ; 
BOUNDS  2.00001  <  XJOIN  <  2.99999; 

MODEL  hormone  =  YJOIN  +  AVESLOPE*(hour-XJOIN)  +  DELSLOPE*ABS(hour-XJOIN); 
DER.YJ0IN=1 ;  DER.aveslope=hour-xjoin;  DER.delslope=abs(hour-xjoin); 
if  hour  >  xjoin   then  DER.xjoin=-aveslope-delslope; 
if  hour  <  xjoin   then  DER.xjoin=-aveslope+delslope; 

OUTPUT  OUT=parsavel  PARMS= YJOIN  AVESLOPE  DELSLOPE  XJOIN  SSE  =  SSE; 
DATA  newl;  SET  parsavel;  FILE  'aiparamsl'; 

SL0PE1=  AVESLOPE-DELSLOPE;    SL0PE2=  AVESLOPE+DELSLOPE; 

IF  KLICK  EQ  1    THEN  PUT  COW  $  (XJOIN  YJOIN)  (11.6)  (SLOPEl  SL0PE2)  (10.5)  SSE  13.8; 

/*     ****  Solve  for  second  interval    ****  */ 
PROC  NLIN  DATA=raw;  BY  cow; 

PARMS  YJOIN=  5  AVESLOPE=  0.1  DELSL0PE=-2      XJOIN=  3.5  ; 
BOUNDS  3.00001  <  XJOIN  <  3.99999; 

MODEL  hormone  =  YJOIN  +  AVESLOPE*(hour-XJOIN)  +  DELSLOPE*ABS(hour-XJOIN); 
DER.YJ0IN=1 ;  DER.aveslope=hour-xjoin;  DER.delslope=abs(hour-xjoin); 
if  hour  >  xjoin   then  DER.xjoin=-aveslope-delslope; 
if  hour  <  xjoin   then  DER.xjoin=-aveslope+delslope; 

OUTPUT  0UT=parsave2  PARMS=YJOIN  AVESLOPE  DELSLOPE  XJOIN  SSE  =  SSE; 
DATA  new2;  SET  parsave2;  FILE  'a:params2'; 

SL0PE1=  AVESLOPE-DELSLOPE;    SL0PE2=  AVESLOPE+DELSLOPE; 

IF  KLICK  EQ  1  THEN  PUT  COW  $  (XJOIN  YJOIN)  (11.6)  (SLOPEl  SL0PE2)  (10.5)  SSE  13.8; 
/*     ****  Solve  for  third  interval    ****      */       Similar  to  above 


XINT  XINT 

Fig.  2A  Fig  2B 

Discussion 

There  are  two  scenarios  regarding  the  general  purpose  of  two-phase  regressions.  In  the  first  a 
researcher  may  have  many,  say  20,  data  sets  from  which  he  wants  to  generate  20  join  points  and  possibly 
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20  pairs  of  slope  parameters  to  be  used  as  data  in  subsequent  analyses.  At  this  point  he  may  wish  to  test 
whether  his  data  fit  the  two-phase  model.  The  obvious  hypothesis  is,  "Does  the  two-phase  model  effect 
a  significant  reduction  in  the  error  sum  of  squares  over  the  simple  linear  model?"  Unfortunately,  the  linear 
model  is  only  one  of  many  possible  subsets  of  the  two-phase  model;  the  quadratic  model  is  another.  A 
test  of  significance  of  the  change  between  the  two-phase  model  and  any  of  its  subsets  is  a  necessary 
but  not  a  sufficient  condition  to  warrant  use  of  the  two-phase  model.  One  may  use  such  a  test  (Overall 
and  Klett,  1972,  p.  422)  as  a  guide,  but  in  practice  any  set  of  data  that  has  even  a  slight  bend  in  the 
regression  line  will  show  a  highly  significant  change  in  favor  of  the  two-phase  model  over  the  linear 
model.  Peiiiaps  the  best  test  is  to  closely  "eye"  the  data's  scattergram  to  see  whether  it  looks  two-phased 
rather  than  linear,  quadratic  or  exponential. 

In  the  second  case  a  researcher  has  combined  his  data  into  one  set  and  after  deriving  the  regression 
parameters  he  may  want  to  test  whether  either  slope  is  different  from  zero.  Nickerson  et  al.  (1989)  pp. 
872-7  describe  such  a  procedure.  Further  study  may  also  be  warranted  for  those  regressions  which 
produce  a  Q.S  solution.  Rather  than  forcing  XJOIN  to  equal  some  Xj,  XJOEN  may  best  be  described  as 
being  the  range  Xj  to  Xj+j  rather  than  a  point.  The  solution  to  this  model  is  given  by  Yeager  and  Ultsch 
(1989). 

Conclusion 

Unconstrained  nonlinear  regression  will  not  in  general  provide  the  correct  solution  a  two-phase 
regression  problem.  Nonlinear  regression  with  bounds  on  the  join  point  values  provides  accurate  answers 
and  does  so  much  more  quickly  and  easily  than  does  the  conventional  approach. 
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Abstract 

The  Multivariate  Analysis  of  Variance  Procedure  (MANOVA)  is  useful  when  a  large  number  of 
response  variables  are  to  be  simultaneously  analyzed  and  joint  effects  are  to  be  determined.  The  test 
criteria  that  correspond  to  the  F-test  in  a  univariate  analysis  are:  Wilks',  Pillais,  Hotelling-Lawley  Trace 
and  Roy's  Greatest  Root.  The  Wilks'  is  an  exact  F-test  and  most  often  used  as  a  test  criterion.  The 
simultaneous  response  of  the  experimental  units  to  all  variables,  considered  as  a  single  response,  generally 
contains  more  information  about  the  total  effect  of  the  treatment  than  does  the  series  of  responses 
considered  individually.  It  is  for  this  reason  that  a  MANOVA  should  precede  a  univariate  test  under 
certain  situations;  an  example  is  presented  and  discussed. 

Introduction 

The  analysis  of  variance  (ANOVA)  is  a  widely  used  statistical  procedure  which  partitions  the 
variation  into  recognizable  sources  and  tests  for  differences  between  the  mean  responses  for  the  treatments. 
In  the  agricultural  sciences,  and  animal  science  in  particular,  it  is  used  in  a  Univariate  Type  of  analysis 
where  one  or  more  main  effects  (X's)  are  tested  in  relation  to  a  single  response  variable  (Y).  This  is  done 
usually  by  comparing  the  variance  of  the  observed  treatment  means  with  an  estimate  of  their  expected 
variance,  if  there  were  no  treatment  differences.  The  response  measure  is  assumed  to  be  normally 
distributed  with  a  constant  variance. 

The  objective  of  this  paper  is  to  introduce  the  Multivariate  ANOVA  using  the  more  widely  used 
Univariate  ANOVA  as  a  basis  and  discuss  briefly  the  application  of  the  Multivariate  ANOVA  in 
agricultural  experiments. 

Multivariate  ANOVA 

Although  the  Univariate  ANOVA  is  adequate  for  studying  one  dependent  variable  such  as  body 
weight,  gain  or  grain  yield  at  a  time,  there  are  situations  when  the  dependents  should  be  studied 
simultaneously.  The  SAS  procedure  has  two  methods  to  calculate  ANOVA 's,  PROC  ANOVA  and  PROC 
GLM.  The  main  advantage  in  using  PROC  ANOVA  is  that  it  is  much  faster,  but  can  only  be  used  for 
balanced  data.  In  the  event  that  the  data  are  unbalanced  and  least  squares  means  are  to  be  derived,  PROC 
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GLM  is  the  procedure  to  use.  In  addition,  PROC  GLM  can  be  used  for  testing  contrasts  (weighted  sums 
of  the  means)  and  will  output  a  data  set  that  contains  residuals.  Unfortunately  GLM  will  not  do  a  means 
separation  such  as  Duncans,  SNK  or  Tukeys  on  least  squares  estimates.  The  same  procedures  (ANOVA 
and  GLM)  can  be  used  for  Multivariate  ANOVAS,  referred  to  as  a  MANOVA. 

The  MANOVA  is  appropriate  under  many  situations.  Oftentimes  there  are  a  large  number  of 
response  variables  and  it  becomes  difficult  to  look  at  each  variable  separately.  Furthermore,  when  there 
is  correlation  between  response  variables  one  needs  to  identify  the  most  important  responses  relative  to 
the  main  effects.  In  instances  where  there  is  a  high  correlation  between  two  responses,  it  could  be  argued 
that  both  responses  are  equally  important  as  shown  through  a  Univariate  analysis  whereas  only  one  may 
be  important.  It  is  for  this  reason  that  a  Multivariate  ANOVA  is  considered  necessary  to  precede  a 
Univariate  analysis  (LitteU  et  al.  1991),  and  only  significant  responses  analyzed  further  by  a  Univariate 
approach.  The  simultaneous  response  of  the  experimental  units  to  all  variables,  considered  as  a  single 
response,  generally  contains  more  information  about  the  total  effect  of  the  treatment  than  does  the  series 
of  responses  considered  separately  (Winer  1971). 

Multivariate  Tests 

In  a  Univariate  ANOVA,  a  main  effect  is  usually  tested  by  a  F  Statistic.  This  is  derived  by  taking 
a  ratio  of  the  sums  of  squares,  the  one  derived  from  the  hypothesis  being  tested  (Treatment  Numerator) 
and  the  other  derived  from  the  unexplained  or  appropriate  error  term  (Denominator).  As  Multivariate 
ANOVAs  have  more  than  one  response  variable  (Yl  to  Yn),  the  sums  of  squares  are  replaced  by  matrices 
of  sums  of  squares  and  cross  products. 

The  Multivariate  linear  model  is  written  as 

Y  =  XB  +  e 

where  Y  is  a  n  x  k  matrix  for  k  dependent/response  variables,  X  is  a  n  x  m  matrix  of  n  observations  on 
M  dependents,  B  is  a  m  x  k  matrix  of  regression  coefficients  and  e  is  the  n  x  k  matrix  of  n  random  errors. 

The  SAS  System  has  four  functions  which  are  used  to  determine  significance  synonymous  with 
the  F-test  in  a  Univariate  analysis.  They  are  Hotelling-Lawley  Trace,  Wilks'  Lambda  criterion,  Pillais 
Trace  and  Roy's  Maximum  Root.  Simulation  criteria  have  not  been  able  to  identify  any  one  of  the  criteria 
as  being  superior  (Littell  et  al  1991)  although  the  Wilks  Lambda  criterion  is  used  by  many  statisticians. 

Example 

Suppose  we  need  to  evaluate  four  swine  feeds:  a  control  (usually  mixed  on  farm),  new  (mixed 
on  farm)  and  Commercial  A  and  Commercial  B.  The  two  dependent  variables  are:  1.  How  long  does 
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it  take  the  pig  to  consume  (TIME)?;  and  2.  How  much  is  eaten  (QUANT)?  The  experiment:  Twenty 
pigs  of  the  same  breed,  sex,  age  and  stage  of  growth  were  randomly  assigned  to  the  four  treatments  C, 
N,  A  &  B  (C  =  Control,  N  =  New,  A  =  Comm  A,  B  =  Comm  B). 
SAS  Code 

DATA  MANOVA; 

*  ONE  WAY  (Rl)  MANOVA  -  SWINE  FEED  TRIAL; 
INPUT  PIGNO  1-3  FEED  $  6-9  TIME  13-14 
QUANT  19-22  .1; 
CARDS; 

PROC  PRINT;  TITLE  'PIG  FEEDING  TRIAL'; 
PROG  GLM; 
CLASS  FEED; 

MODEL  TIME  QUANT  =  FEED; 

MEANS  FEED; 

MANOVA  H  =  FEED; 

OUTPUT  OUT  =  RESIDS  R  =  Rl  R2; 

RUN; 

Output 

The  order  in  which  these  tables  are  presented  is  not  necessarily  the  order  in  which  the  SAS  output 
is  obtained.  This  has  been  done  to  show  you  in  what  order  the  output  should  be  read. 

Table  1.         Eigenvalues  and  Vectors 


General  Linear  Models  Procedure 
Multivariate  Analysis  of  Variance 


Characteristic  Roots  and  Vectors  of:  E  Inverse  *  H,  where 
H  =  Type  III  SS&CP  Matrix  for  FEED 


E  =  Error  SS&CP  Matrix 


Characteristic 
Root 


Percent 


Characteristic  Vector 
TIME  QUANT 


V'EV=1 


2.3315933541 
0.0904626730 


96.27 
3.73 


0.01453729 
0.07588082 


0.50402729 
0.04497715 
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The  characteristic  roots  are  the  eigenvalues.  These  represent  linear  functions  of  the  dependent  or 
response  variables.  The  percent  96.27  represents  the  variation  that  each  linear  function  explains.  The 
vectors  with  large  values  show  which  of  the  two  response  variables  are  explained  more  by  the  model. 
This  is  sometimes  very  useful.  In  our  example,  QUANT  =  0.504  is  explained  more  than  TIME  =  0.014. 

Table  2  shows  the  MANOVA  analysis  testing  the  null  hypotheses  of  NO  overall  FEED  effect 

Table  2  MANOVA  analysis  for  FEED  with  test  criteria 

H  =  Type  III  SS&CP  Matrix  for  FEED  E  =  Error  SS&CP  Matrix 


S=2 

M=0 

N=6.5 

Statistic 

Value 

F 

Num  DF 

Den  DF 

PR>F 

Wilks'  Lambda 

0.275256 

4.5302 

6 

30 

0.0022 

Pillai's  Trace 

0.782801 

3.43 

6 

32 

0.0100 

Hotelling-Lawley  Trace 

2.422056 

5.6515 

6 

28 

0.0006 

Roy's  Greatest  Root 

2.331593 

12.435 

3 

16 

0.0002 

Note:    F  Statistic  for  Roy's  Greatest  Root  is  an  upper  bound. 
Note:    F  Statistic  for  Wilks'  Lambda  is  exact. 

There  are  four  Multivariate  tests.  Using  Wilks',  which  is  an  exact  F-test  with  6  degrees  of 
freedom  for  the  numerator  and  30  for  the  denominator  our  probability  of  rejecting  is  0.0022.  The 
conclusion  is  that  FEED  affects  TIME  and  QUANT  jointly.  The  same  output  may  be  obtained  by  using 
contrast  statements  as  shown  below: 

FEED  100-1, 
FEED  10-10, 
FEED  1-10  0; 

The  probability  for  the  hypothesis  of  no  overall  FEED  effect  (HJ  can  sometimes  be  smaller  for  joint 
effects  compared  to  independent  (Univariate)  effects.  Therefore,  one  could  detect  significant  differences 
by  viewing  variables  together  that  you  would  not  detect  by  looking  at  variables  individually.  This  is  an 
advantage  in  the  MANOVA  approach. 

Table  3  shows  the  mean  and  descriptive  statistics  generated  due  to  inclusion  of  the  PROC  MEANS 
statement  in  the  code. 
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Table  3 


Descriptive  Statistics 


General  Linear  Models  Procedure 


Level  of   TIME  

FEED  N                    Mean  SD 

ComA  5                  26.2000000  4.14728827 

ComB  5                  27.2000000  3.42052628 

Com  5                  25.2000000  3.11448230 

New  5                  27.8000000  2.58843582 

Level  of   QUANT  

FEED  N                    Mean  SD 

ComA  5                   7.96000000  0.36469165 

ComB  5                   7.70000000  0.53851648 

Com  5                   9.44000000  0.46151923 

New  5                   8.80000000  0.64807407 


The  data  shows  greater  differences  in  QUANT  compared  with  TIME.  This  will  become  evident 
with  the  Univariate  tests. 

Tables  4  and  5  show  the  Univariate  tests  for  TIME  and  QUANT. 


Table  4 


Univariate  test  for  TIME  with  FEED  as  Main  Effect 


Dependent  Variable:  TIME 


Source 
Model 
Error 

Corrected  Total 


DF 

3 

16 
19 


Sum  of  Squares 
19.60000000 
181.20000000 
200.80000000 


F  Value 
0.58 


PR>F 
0.6385 


R-Square 
0.097610 


C.V. 

12.65137 


TIME  Mean 
26.60000000 


Source 
FEED 


DF 

3 


Type  I  SS 
19.60000000 


F  Value 
0.58 


PR  >  F 
0.6385 


Source 
FEED 


DF 

3 


Type  III  SS 
19.60000000 


F  Value 
0.58 


PR  >  F 
0.6385 
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Table  5 


Univariate  tests  for  QUANT  with  FEED  as  Main  Effect 


General  Linear  Models  Procedure 


Dependent  Variable:  QUANT 


Source 
Model 
Error 

Corrected  Total 


DF 

3 

16 
19 


Sum  of  Squares 
9.51350000 
4.22400000 
13.73750000 


F  Value 
12.01 


PR>F 
0.0002 


R-Square 
0.692520 


C.V. 

6.062647 


TIME  Mean 
8.47500000 


Source 
FEED 


DF 

3 


Type  I  SS 
9.51350000 


F  Value 
12.01 


PR>F 
0.0002 


Source 
FEED 


DF 

3 


Type  III  SS 
9.51350000 


F  Value 
12.01 


PR  >  F 
0.0002 


*  Note:  SAS  gives  you  type  I  and  III  SS  by  default. 


FEED  is  not  a  significant  source  of  variation  for  TIME,  P  =  0.6385,  but  FEED  is  highly 
significant  for  QUANT.  This  same  idea  was  noted  by  looking  at  the  eigenvalues  in  Table  1.  However, 
although  QUANT  is  affected  more  by  FEED,  both  TIME  AND  QUANT  are  affected  joindy  by  FEED 
(Table  2). 

There  are  option  statements  in  the  procedure  that  are  helpful:  the  PRINTH  and  PRINTE,  print  the 
hypothesis  and  error  matrices,  respectively.  The  latter  provides  partial  correlations  from  the  sums  of 
squares  and  cross  products  that  are  adjusted  for  all  model  effects.  The  SHORT  option  gives  a  condensed 
printout.  If  the  data  is  balanced,  PROC  ANOVA  may  be  used  but  PROC  GLM  should  be  used  for 
unbalanced  data.  If  CONTRAST  statements  are  to  be  used,  one  should  use  PROC  GLM.  A  factorial 
MANOVA  will  be  discussed  at  the  Workshop. 

The  MANOVA  much  like  the  ANOVA  assumes  that  the  response  variables  and  residuals  are 
normally  distributed.  The  following  code  can  be  linked  to  the  code  provided  on  page  2,  as  a  new 
OUTPUT  data  set  has  been  created  for  variables  Rl  (TIME)  and  R2  (QUANT).  The  probabilities  of 
rejecting  a  null  hypothesis  was  0.134  for  Rl  and  0.07  for  R2  suggesting  that  both  response  variables  did 
not  deviate  significantly  from  normality. 
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SAS  Code 

PROC  UNIVARIATE  Plot  Nomial  data  =  RESIDS; 
VAR  Rl  R2; 

Title  'Check  normality  assumptions'; 
RUN; 

Multivariate  tests,  although  they  lack  power,  are  useful  to  detect  joint  effects.  PROC  ANOVA 
and  PROC  GLM  can  be  used  to  perform  MANOVAS.  Although  traditional  approaches  have  looked  at 
Univariate  tests  only,  it  may  be  better  in  certain  in  instances  to  look  at  MANOVA  output  and  then  decide 
on  what  Univariate  variables  are  significant. 
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The  Determination  of  Risk  of  Disease  in  Cattle  Feedlots: 
A  Case-Control  Study 

Casey  Schipper 
Animal  Health  Division 
Alberta  Agriculture 
Edmonton,  Alberta 

Introduction 

Shipping  fever  is  a  major  cause  of  sickness  in  feedlot  calves.  The  highest  incidence  of  this  disease 
occurs  within  six  weeks  after  arrival  at  the  lot.  A  variety  of  infectious  agents  as  well  as  a  multitude  of 
management  factors  are  now  believed  to  be  responsible  for  causing  this  disease.  Factors  such  as  mixing 
cattle  from  different  sources,  early  silage  feeding,  processing,  vaccination  and  water  medication  were 
associated  with  increased  deaths  and  treatment  costs  in  Ontario  feedlots  (1). 

Objectives 

The  objective  of  this  workshop  presentation  is  to  explain  how  simple  cross-tabulation  and  logistic 
regression  methods  may  be  used  to  estimate  the  relative  importance  of  contributing  management  factors 
(exposure  factors)  on  the  risk  (probability  of  occurrence)  of  a  given  outcome  variable  (shipping  fever). 

Materials  and  Methods 

Each  line  of  this  fictitious  data  set  represents  a  partial  health  history  of  each  of  2660  autunm 
weaned  steer  calves.  The  case-control  study  consisted  of  284  calves  affected  with  shipping  fever  and  2376 
unmatched  controls.  The  calves'  individual  eartag  numbers  were  drawn  by  simple  random  selection  from 
the  records  of  a  cattle  feedlot  in  Alberta,  housing  approx.  6000  calves.  The  calves  (units  of  analysis)  in 
this  sample  population  weighed  between  180  and  340  kg  (400  and  750  lbs),  and  entered  the  feedlot 
between  Nov.  15  and  Dec  20,  1990. 

For  the  purpose  of  this  analysis  the  following  exposure  variables  were  considered: 

1.  SRC  -  About  10%  of  the  sampled  calves  had  been  purchased  through  auction 
markets  in  Manitoba. 

2.  DH  -  About  36%  of  all  calves  required  dehorning  upon  arrival. 

3.  DC  -   About  22%  of  the  sample  population  underwent  castration  as  well  as 
dehorning  at  that  time. 

4.  VAC  -  Vaccination  against  shipping  fever  was  given  to  36%  of  the  calves. 

5.  SUR  -  Processing  included  either  castration  or  dehorning  or  both. 

6.  BWT  -  Midpoint  body  weights  selected  were  438,  513,  588,  663  and  738  lbs. 
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The  outcome  variable  "Shipping  fever"  was  defined  as  a  condition  diagnosed  by  the  feedlot 
penchecker/treatment  person.  Criteria  for  diagnosis  were: 

-  body  temperature  greater  than  40.5°C  or  (105°F) 

-  off  feed,  depression,  abnormal  breathing 

Statistical  Analysis 

The  strength  of  association  between  an  exposure  factor  and  disease  is  said  to  be  the  relative  risk 
(RR),  calculated  as  the  ratio  of  the  rate  of  disease  in  the  exposed  group  to  that  in  the  unexposed  group. 
The  RR  =  1  if  there  is  no  association  between  factor  and  disease.  The  further  the  RR  deviates  from  1  the 
greater  is  the  strength  of  association. 

The  odds  ratio  (OR)  is  used  in  the  analysis  of  case-control  studies,  since  it  closely  approximates 
the  value  of  the  RR  estimate.  The  OR  is  derived  from  the  ratio  a*d/b*c,  where, 
a  =  number  of  cases  of  disease  in  the  exposed  group 
b  =  number  of  exposed  controls 
c  =  number  of  unexposed  cases  of  disease 
d  =  number  of  unexposed  controls 
The  calculation  of  approximate  95%  confidence  limits  on  the  OR  is  based  on  the  natural  log  Qn) 
transformation  of  the  approximate  limits  of  In(OR).  One  needs  to  calculate  the  variance  of  the  OR,  as 
foUows  (2) 

var  (InOR)  =  (1/a  +  1/b  +  1/c  +  1/d) 
Since  In(OR)  has  a  normal  distribution  in  large  samples  the  approximate  95%  confidence  limits 
for  hi(OR)  are, 

hi(OR)  +  1.96  *  (1/a  +  1/b  +  1/c  +  1/df  ^ 
The  approx.  95%  lower  and  upper  limits  of  In(OR)  are  the  antilog  of  the  lower  and  upper  limits 
of  In(OR)  ie. 

ORl  =  OR  exp[-1.96  *  {var  (lnOR)}f  ^]  and 

ORu=  OR  exp[+1.96  *  {var  (lnOR)}f  ^] 
One  needs  to  be  cautious  in  interpreting  risk  of  disease  associated  with  exposure  factors.  The 
apparent  or  observed  association  may  actually  be  due  to  another  variable.  The  term  confounding  refers 
to  the  effect  of  an  extraneous  variable  that  wholly  or  partially  accounts  for  the  apparent  effect  of  the  study 
exposure  factor.  A  confounding  variable  by  definition  is  associated  with  the  main  independent  variable 
studied  and  the  outcome  (disease)  variable,  but  is  not  a  consequence  of  exposure  (3).  Confounding 
variables  can  also  be  determinants  of  disease.  If  confounders  are  ignored  in  the  analysis  the  real 
association  between  the  exposure  factor  of  interest  and  disease  may  be  distorted.  The  control  of 
confounding  variables  may  be  obtained  by  analyfic  methods  (4). 
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The  statistical  method  used  for  summarizing  associations  in  multiple  tables,  is  known  as  the 
Mantel-Haenzel  (MH)  technique.  Different  2-by-2  tables  showing  number  of  cases  and  controls  for  each 
level  of  an  exposure  factor  are  constructed  for  each  value  level  of  the  confounding  variable.  A  limitation 
of  the  MH  technique  is  tiie  sheer  numbers  of  2-by-2  tables  that  need  to  be  constructed  to  account  for  each 
category  of  perhaps  many  confounding  variables  (5).  For  instance  in  tiiis  study  2*  =  256  tables  would 
have  been  required  to  obtain  a  Summary  Odds  Ratio  for  the  risk  of  shipping  fever. 

The  data  set  was  created  in  Lotus  1-2-3.  The  Lotus  data  file  was  imported  and  converted  into  a 
epidemiologic  data  analysis  program  "EPI5".  The  ANALYSIS  component  of  this  program  is  designed 
to  do  frequencies,  cross  tabulations,  means,  graphs,  and  linear  regression.  The  program  computes  OR 
values  and  their  95%  confidence  intervals  and  enables  one  to  calculate  Mantel  and  Haenzel  Summary 
Odds  ratios.  The  STATCALC  module  calculates  statistics  from  tables  entered  from  the  keyboard.  It  also 
performs  calculations  for  single  and  stratified  2-by-2  tables,  determines  sample  size  for  cross-sectional, 
cohort  and  case-control  studies  and  provides  single  and  stratified  trend  analysis  (6). 

When  the  potential  effects  of  a  large  number  of  suspected  concurrent  exposure  factors  need  to  be 
evaluated,  multivariate  logistic  regression  is  commonly  used.  Multivariate  modeUing  of  relative  risk 
functions  has  a  20-year  history  in  epidemiologic  research.  The  logistic  regression  approach  is  generally 
used  in  case-control  studies  of  chronic  disease.  It  is  also  useful  in  the  analysis  of  disease  risk  in  catUe 
feedlots  (7).  The  technique  allows  one  to  examine  the  effect  of  one  exposure  or  management  factor  while 
the  values  of  other  variables  in  the  regression  equation  are  held  constant  mathematically. 

MULTLR  is  a  Pascal  microcomputer  program,  which  performs  interactive,  menu-driven  multiple 
logistic  regression  analyzed  using  conditional  and  unconditional  maximum  likelihood  estimation  methods. 
The  program  calculates  most  standard  statistics  and  allows  factoring  of  categorical  or  continuous  variables 
by  two  distinct  methods  of  contrast.  A  built-in,  descriptive  statistics  option  allows  the  user  to  inspect  the 
distribution  of  cases  and  controls  across  categories  of  any  given  variable  (8). 

The  binary  outcome  variable  D  (disease)  is  conditionally  related  to  a  series  of  independent 
regressors  (discrete  or  continuous  exposure  variables)  Xj,  X2,  ...,Xp  through  the  formula  (9): 

Prob(D  =  II  Xj,  X2,     Xp)  =  {1  +  exp[  -(a  +  bi*Xi  +  b2*X2+  ...  +  bp*Xp)]}. 
where  D  denotes  the  presence  (D  =  1)  or  absence  (D  =  0)  of  disease,  and  x  represents  a  set  of  p  variables 
X  =  (Xi,  X2,  X3,...Xp).   The  X  variables  represent  any  potential  risk  factor,  confounding  factor  or  any 
interaction  term  of  interest.  The  b  coefficients  are  parameters  that  represent  the  effects  of  the  x's  on  the 
risk  (probabihty)  of  disease. 

The  coefficient  for  any  variable  included  in  the  logistic  regression  model  depends  on  the  entire 
set  of  variables  included  in  the  model.  Thus  a  variable  that  by  itself  may  be  significantiy  associated  with 
disease,  may  have  a  coefficient  that  renders  it  insignificant  in  multivariate  analysis. 
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One  should  ensure  that  in  judging  the  suitability  of  a  fitted  model  the  estimated  effects  of  the 
variables  included  in  the  analysis  make  biological  sense.  Logistics  regression  models  are  used  to  obtain 
point  (mean)  estimates  and  suitable  (95%)  confidence  limits  for  the  probability  (risk)  of  disease  and 
unconditional  or  conditional  odds  ratios. 

A  detailed  discussion  on  the  rationale,  methodology  and  interpretation  of  logistic  regression 
functions  may  be  found  in  suitable  textbooks  on  the  subject  (10). 

Results 

The  exposure  factor  with  the  greatest  unconditional  impact  on  the  risk  of  shipping  fever  was  the 
application  of  a  surgical  procedure,  ie.  either  dehorning  and  castration  (OR  =  4.63,  95%  =  3.55,  6.04). 
The  source  or  origin  of  the  cattle  also  had  a  significant  impact  on  disease  unconditionally  (OR  =  1.96, 
95%  CL  =  1.33,  2.89).  The  weight  of  calves,  categorized  in  5  groups,  was  significantly  associated  with 
the  risk  of  disease.  Animals  weighing  between  400  and  475  lbs  had  a  disease  risk  of  19.2%  while  only 
1.8%  of  those  weighing  more  than  700  lbs  needed  treatment  for  shipping  fever.  Dehorning  by  itself  was 
not  associated  with  an  increased  risk  of  illness  (p  =  0.71).  Similarly,  vaccination  against  Bovine 
Respiratory  disease  had  no  sparing  effect  (p  =  0.58)  on  the  risk  of  shipping  fever  in  this  data  set  (Table 
1). 

With  respect  to  castration  and  dehorning,  the  source  of  calves  was  a  confounder  (p  =  0.002),  but 
vaccination  and  body  weight  were  not.  With  respect  to  source  however,  body  weight  was  a  significant 
confounding  variable  (X^  =  1 10.2  with  4  degrees  of  freedom). 

The  full  logistic  regression  model  only  included  dehorning/castration  procedure,  body  weight,  and 
the  interaction  term  source*dehoming/castration  as  significant  disease  risk  determinants  (p  <  0.05).  The 
effect  of  source  became  insignificant  when  the  remaining  factors  in  the  equation  were  controlled 
statistically.  There  was  no  interaction  between  surgery  and  the  weight  of  calves  (Table  2). 

Taken  in  combination,  the  major  factors  influencing  the  probability  of  disease  were  light  weight 
calves,  dehorning/castration  and  the  source  of  calves  (Table  3). 

Conclusion 

The  best  calves  to  purchase  from  a  disease  control  point  of  view  based  on  this  sample  data  set, 
are  dehorned  Alberta  yearling  steers.  More  than  one  half  (64%)  of  all  light  weight  Manitoba  homed  bull 
calves,  entered  into  an  Alberta  feedlot  under  conditions  prevalent  in  this  data  set  would  likely  require 
treatment  for  shipping  fever. 
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Table  1.  Calf  Sample  Distribution  of  Shipping  Fever  Risk  in  Fictitious  Alberta  Feedlot 

(November  15,  1990  to  January  30,  1991) 


Variahlf  TF) 

Variable 

Tntpmrptatinn 

Number  (%) 

P'YDOQPii 

Number  (%) 

Source 

0 

Alberta 

2347  (91.6) 

244(10.0) 

1 

Manitoba 

223  (8.4) 

40  (17.9) 

Dehorning 

0 

no  dehorning 

1706  (64.1) 

185  (10.8) 

1 

dehorning 

954 

99  (10.4) 

Dehorning  & 
Castration 

0 

No 

2084  (78.3) 

140  (6.7) 

1 

1 

576  (1\  1^ 

144  Vi\ 

V  accinaiion 

U 

INO 

1 

I  es 

QQ  /I  n  o\ 

Body  Weight 

438 

lbs 

950  (35.7) 

182  (19.2) 

513 

lbs 

767  (28.8) 

68  (8.9) 

588 

lbs 

605  (22.7) 

22  (3.6) 

663 

lbs 

224  (8.4) 

10  (4.5) 

738 

lbs 

114  (4.3) 

2  (1.8) 
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Table  2.  Maximum  Likelihood  Estimates  of  Logistic  Parameters  Relating  Three  Risk  Factors 
to  the  Development  of  Shipping  Fever  in  the  Feedlot  (November  13,  1990  to 
January  30,  1991). 


Parameter       Estimate  b       Std.  Error 


Variable 

(b) 

coeff 

of  coeff 

Odds  Ratio 

P-Value 

Xq  Intercept 

bo 

-2.252 

0.095 

- 

0.000 

Xi  Source 

bi 

0.226 

0.365 

1.25 

0.536 

X2  Castrate  &  dehorn 

b2 

1.346 

0.149 

3.84 

0.000 

X3  Body  Weight 

b3 

-1.899 

0.307 

0.15 

0.000 

Xi  *  X2 

C^ 

1.256 

0.473 

3.51 

0.008 

X2  *  X3 

C2 

0.275 

0.416 

1.31 

0.509 

Maximized  log  Likelihood  L  (b)  =  -779.38 

Estimated  Covariance  Matrix 

bo 

b, 

bz 

b3 

C2 

bo  0.009 

bi  -0.007 

0.133 

b2  -0.009 

0.007 

0.022 

b3  -0.008 

-0.010 

0.008 

0.094 

Ci  0.008 

-0.133 

-0.016 

0.011 

0.224 

C2  0.008 

0.011 

-0.018 

-0.094 

-0.044 

0.1773 
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Table  3.  A  Summary  of  Maximum  Likelihood  Predicted  Effects  of  Significant  Exposure 
Factors  and  the  Risk  of  Shipping  Fever  in  an  Alberta  Feedlot  (November  15,  1990 
to  January  30, 1991). 


Outcome  Variable 

Values  of  Exposure 
Variables^ 

Predicted  Treatment 
Rates  Point  Estimate 

(%) 

Predicted  Odds  Ratios 
Point  Estimate 

Shipping  Fever 

SRC^ 

DC^ 

WT* 

0 

0 

1 

1.6 

0.15  (6.67) 

1 

0 

1 

1.9 

0.19  (5.26) 

0 

1 

1 

7.4 

0.75  (1.33) 

0 

0 

0 

9.5 

1.00  (1.00) 

1 

0 

0 

11.6 

1.25 

1 

1 

1 

25.9 

3.33 

0 

1 

0 

28.8 

3.84 

1 

1 

0 

64.0 

16.91 

^Values  adjusted  for  the  SRC*DC  and  DC*WT  interactions. 

=  Alberta  cattie  1  =  Manitoba  cattle 

^0  =  no  surgical  stress  1  =  dehorned  and  castrated  on  arrival 

"^0  =  calves  weighing  less  than  550  lbs  1  =  calves  weighing  550  lbs  and  more 
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Abstract 

Statistical  analysis  plays  a  vital  role  in  environmental  restoration  (ER)  activities.  This  paper  is 
an  attempt  to  show  where  statistics  fits  into  the  ER  process.  The  statistician,  as  member  of  the  ER 
planning  team,  works  collaboratively  witii  the  team  to  develop  tiie  site  characterization  sampling  design, 
so  that  data  of  the  quality  and  quantity  required  by  the  specified  data  quality  objectives  (DQOs)  are 
obtained.  At  Uie  same  time,  the  statistician  works  with  the  rest  of  the  planning  team  to  design  and 
implement,  when  appropriate,  the  observational  approach  to  streamline  the  ER  process  and  reduce  costs. 
The  statistician  will  also  provide  Uie  expertise  needed  to  select  or  develop  appropriate  tools  for  statistical 
analysis  that  are  suited  for  problems  that  are  common  to  waste-site  data.  These  data  problems  include 
highly  heterogeneous  waste  forms,  large  variability  in  concentrations  over  space,  correlated  data,  data  that 
do  not  have  a  normal  (Gaussian)  distribution,  and  measurements  below  detection  limits.  Other  problems 
include  environmental  transport  and  risk  models  that  yield  highly  uncertain  predictions,  and  the  need  to 
effectively  communicate  to  the  public  highly  technical  information,  such  as  sampling  plans,  site 
characterization  data,  statistical  analysis  results,  and  risk  estimates.  Even  though  some  statistical  analysis 
methods  are  available  "off  the  shelf  for  use  in  ER,  these  problems  require  the  development  of  additional 
statistical  tools,  as  discussed  in  this  paper. 

Introduction 

The  U.S.  Department  of  Energy  (DOE)  has  estimated  that  it  will  cost  many  billions  of  dollars  to 
clean  up  its  large  volume  of  hazardous,  radioactive,  and  mixed  waste  (DOE,  1990).  As  a  consequence, 
the  DOE  is  taking  a  close  look  at  ways  to  reduce  environmental  restoration  (ER)  costs  while  still  ensuring 
the  scientific  validity  of  the  ER  process.  Achieving  scientific  validity  begins  by  developing  thorough  plans 
for  characterizing  the  waste  site.  These  plans  must  specify  both  the  quality  and  the  quantity  of  site- 
characterization  data  required  .  That  is,  it  is  necessary  to  specify  how  much  uncertainty  in  the  nature  and 
extent  of  contamination  at  the  site  can  be  allowed.  Then  efficient  field  sampling  plans  that  will  generate 
sufficient  data  of  required  quality  can  (and  must)  be  developed. 
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In  this  paper,  we  discuss  statistical  issues  and  needs  related  to  the  planning  and  design  of  site 
characterization  sampling  studies  and  to  the  analysis  of  generated  data.  We  begin  by  discussing  how  to 
deal  with  uncertainty  (and  cost)  in  the  site-characterization  process. 

Dealing  with  Uncertainty  in  Site  Assessments 

It  is  impossible  to  remove  all  uncertainty  about  the  nature,  magnitude,  and  spatial  patterns  of 
contamination  at  a  hazardous  waste  site,  no  matter  how  thorough  sampling  and  measurement  for  the  site 
characterization  effort  may  be.  We  must  therefore  deal  with  the  issue  of  uncertainty,  which  in  turn  means 
that  we  must  deal  with  statistical  concepts  and  metiiods  that  are  related  to  the  design  of  sampling  programs 
so  that  we  can  quantify  and/or  reduce  uncertainty.  To  deal  with  uncertainty  requires  careful  planning  even 
before  any  samples  are  collected.  An  important  aspect  of  this  planning  process  is  the  specification  of  data 
quality  objectives  (DQOs). 

These  DQOs  are  qualitative  and  quantitative  statements  tiiat  specify  tiie  quality  of  data  that  must 
be  obtained.  They  are  a  tool  to  answer  the  questions:  "What  type  and  quality  of  data  are  needed  to  answer 
key  questions?"  and  "How  do  we  know  when  we  have  enough  data?"  (Neptune,  et  al.  1990).  Specific 
DQOs  are  determined  by  the  end  use  of  the  data  and  so  are  established  to  ensure  that  the  data  are  both 
sufficient  and  of  adequate  quality  for  their  intended  use.  Once  DQOs  are  established  by  the  ER  planning 
team,  they  can  be  used  by  tiie  statistician  in  developing  a  sampling  design  that  wiU  yield  the  necessary 
data  at  minimum  cost. 

The  steps  in  establishing  DQOs  are  as  follows  (Neptune,  et  al.  1990;  EPA  1987a;  EPA  1987b): 

1 .  Carefully  state  the  problem  to  be  addressed  or  the  decision  to  be  made. 

2.  Identify  the  information  required  to  select  an  appropriate  course  of  action. 

3.  Articulate  the  specific  role  that  data  will  play  in  selecting  the  course  of  action. 

4.  Specify  the  type  of  data  needed. 

5.  Specify  the  way  the  data  will  be  used. 

6.  Specify  (by  means  of  an  iterative  process  that  involves  both  the  decision-maker  and  technical 
support  staff)  the  degree  of  certainty  desired  in  the  conclusions  to  be  derived  from  the  data. 

7.  Optimize  the  sampling  design  for  data  collection  to  achieve  the  required  degree  of  certainty  in  the 
conclusions  at  minimum  cost. 

Neptune,  et  al.  (1990)  have  illustrated  the  procedure  with  a  case  study.  The  question  addressed 
in  the  case  study  is  whether  a  site  at  which  railroad  ties  and  creosote-soaked  timbers  were  stored  and 
burned  posed  an  unacceptable  risk  to  site  workers  and  visitors  as  a  result  of  exposure  to  polyaromatic 
hydrocarbons  (PAH)  in  soil.  After  carefully  considering  the  human  health  consequences  and  consulting 
with  the  toxicologist  and  site  engineers,  the  project  manager  assigned  acceptable  decision  error  rates  for 
various  risk  levels  and  associated  levels  of  average  PAH  soil  concentrations.  Then  a  soil  sampling  plan 
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was  developed  that  was  expected  to  achieve  these  DQOs  at  minimum  cost.  The  sampling  plan  specified 
the  number  of  composite  soil  samples  that  should  be  collected  to  estimate  average  soil  PAH  concentrations 
with  sufficient  precision  to  achieve  the  DQOs,  i.e.,  such  that  the  specified  decision  error  rates  were  not 
exceeded. 

Dealing  with  Costs  Using  the  Observational  Approach 

One  lesson  learned  from  ten  years  of  experience  of  the  U.S.  Environmental  Protection  Agency 
(EPA)  with  conducting  clean-up  actions  under  the  Comprehensive  Environmental  Response,  Compensation 
and  Liability  Act  (CERCLA)  is  that  the  clean-up  process  must  be  streamlined  to  avoid  high  costs  and  long 
delays  in  conducting  remedial  actions  (EPA  1989a).  The  DOE  and  the  EPA  are  currently  evaluating  the 
"Observational  Approach"  as  a  framework  for  streamlining  the  clean-up  process  while  managing  the 
uncertainty  inherent  in  site  assessment.  Overviews  of  the  method  have  been  published  by  Smyth  and 
Quinn  (1991)  and  by  Myers  and  Gianti  (1989).  The  advantages  and  limitations  of  the  method  are 
discussed  by  Peck  (1969).  Brown,  et  al.  (1989)  discuss  the  application  of  the  method  to  the  remediation 
of  hazardous  waste  sites.  The  method  is  a  way  of  initiating  remedial  action  at  a  waste  site  without  full 
characterization  of  the  nature  and  extent  of  contamination.  The  observational  approach  is  intended  to 
reduce  costs  by  accepting  greater  uncertainty  and  allowing  earlier  selection  of  a  remedial  action  approach 
based  on  probable  conditions  at  the  waste  site. 

The  uncertainty  in  knowledge  about  conditions  at  the  site  is  taken  into  account  by  making 
contingency  plans  for  handling  deviations  from  probable  conditions  if  they  occur  during  remedial  action. 
In  essence,  the  observational  approach  requires  conducting  thorough  up-front  planning  to  identify 
uncertainties  and  determine  both  possible  and  probable  conditions  at  the  site.  The  remedial  action  program 
is  designed  for  probable  conditions,  but  contingency  plans  are  prepared  in  case  deviations  from  the 
probable  conditions  occur  during  remedial  action.  As  indicated  by  Smyth  and  Quinn  (1991),  the  DOE 
has  endorsed  the  concept  (DOE,  1990),  and  the  EPA  has  endorsed  an  equivalent  approach  (EPA,  1989b). 

Statistical  Needs  in  Environmental  Restoration 

It  is  generally  recognized  that  statistical  methods  should  be  used  in  ER  projects.  However,  a 
number  of  problems  associated  with  using  statistical  methods  must  be  addressed.  The  following  is  a  short 
discussion  of  some  of  these  problems. 
Understanding  the  Role  of  the  Statistician 

A  statistician  can  contribute  to  the  ER  process  in  many  ways,  including  helping  to  define  the 
DQOs,  developing  sampling  plans,  conducting  data  analyses,  reviewing  analytical  laboratory  protocols  and 
quality  control  procedures,  developing  strategies  for  risk  assessment,  quantifying  the  uncertainty  of  model 
predictions  by  using  uncertainty  and  sensitivity  analyses,  developing  graphics  for  communicating  data  and 
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results,  reviewing  draft  reports,  and  contributing  to  peer-reviewed  papers  and  reports.  These  contributions 
will  be  most  effective  when  the  statistician  is  a  member  of  the  ER  planning  team,  preferably  Irom  the  very 
beginning  of  the  ER  planning  process. 
Developing  Data  Quality  Objectives 

Data  quality  objectives  must  be  specified  as  part  of  the  planning  process.  If  this  is  not  done,  no 
one  will  know  when  to  stop  collecting  data,  and  a  lot  of  unnecessary  or  inappropriate  data  are  likely  to 
be  the  result.  This  idea  of  specifying  a  priori  the  quantity  of  data  that  is  needed  is  a  familiar  concept  to 
statisticians.  For  example,  the  method  statisticians  use  to  determine  the  number  of  samples  required  for 
estimating  a  mean  involves  specifying  the  required  accuracy  of  the  estimated  mean  and  the  confidence 
required  in  achieving  that  accuracy  (Gilbert,  1987,  pp.  30-42).  The  entire  ER  plaiming  team  needs  to 
understand  and  support  the  use  of  DQOs  and  to  take  part  in  determining  what  those  DQOs  should  be. 
The  DQO  approach  will  provide  the  information  needed  to  develop  efficient  sampling  plans  and  associated 
statistical  analysis  methods. 
Reducing  Cost 

Ways  must  be  found  to  reduce  the  cost  of  ER.  Efficient  planning  via  the  Observational  Approach 
and  the  specification  of  DQOs  are  two  important  tools  for  that  purpose.  But  other  methods  are  also 
available.  For  instance,  the  compositing  of  several  field  samples  into  one  thoroughly  mixed  sample,  which 
is  then  subsampled  for  analysis,  can  reduce  analytical  costs  when  the  method  is  applicable  (Gilbert,  1987; 
Bolgiano,  et  al.  1990).  For  example,  concentrations  of  composite  soil  samples  were  used  to  evaluate  the 
need  for  further  removal  of  soil  at  sites  contaminated  with  dioxin  in  the  State  of  Missouri  (Exner,  et  al. 
1985).  The  use  of  in  situ  detectors  in  place  of  a  portion  of  the  environmental  samples  can  also  sometimes 
reduce  costs.  For  example,  in  situ  spectrometry  was  used  in  the  United  States  on  the  Nevada  Test  Site, 
a  nuclear  weapons  testing  area,  to  estimate  the  spatial  distribution  and  total  inventory  of  the  important 
anthropogenic  radionuclides  in  the  surface  soil  (Mc Arthur,  1991).  When  in  situ  detectors  are  used  to 
estimate  environmental  concentrations,  it  is  usually  necessary  to  quantitatively  relate  the  in  situ  detector 
readings  to  concentrations  in  samples  collected  at  the  same  locations.  This  quantification  requires  a 
statistical  analysis  of  both  the  in  situ  and  the  sample  data  collected  using  a  valid  statistical  design 
according  to  established  DQOs. 
Coping  with  Spatial  Variability 

Contaminant  concentrations  can  vary  greatly  over  space.  Site  characterization  usually  implies 
estimating  what  contaminants  are  present,  where  they  are  located,  and  their  concentrations.  If  contaminant 
concentrations  are  not  too  heterogeneous  over  space,  geostatistical  techniques,  such  as  kriging,  can  be  used 
to  estimate  the  spatial  pattern  by  taking  into  account  spatial  correlations  (Flatman,  1984;  Gilbert  and 
Simpson,  1985).  The  use  of  geostatistics  for  evaluating  the  attainment  of  clean-up  standards  is  discussed 
in  EPA  (1989c).  However,  the  presence  of  hot  spots  is  another  complicating  problem.  Although  simple 


78 


methods  are  available  for  determining  the  spacing  between  points  on  a  sampling  grid  required  to  detect 
hot  spots  with  specified  probability  (Gilbert,  1987),  the  number  of  sampling  locations  required  can  be 
prohibitively  large.  Sometimes  it  is  possible  to  reduce  variability  by  compositing  and  mixing  samples. 
Also,  in  situ  detector  measurements  can  be  less  variable  than  those  of  single  small  samples  because  the 
in  situ  measurements  measure  relatively  large  volumes  of  soil.  However,  such  smoothing  out  of  sample 
spatial  variabilities  can  also  hide  hot  spots. 
Characterizing  Highly  Heterogeneous  Waste 

Many  hazardous  waste  sites  contain  heterogeneous  materials;  e.g.,  concrete,  clothing,  liquids, 
bottles,  tires,  paper,  and  shredded  autos  may  all  occur  at  a  site.  These  materials  may  be  packed  in  barrels 
or  lying  loose.  A  given  barrel  may  contain  many  different  types  of  waste,  unknown  unless  the  barrel  is 
opened  and  inspected,  which  is  an  expensive  and  possibly  dangerous  operation.  The  EPA  is  searching 
for  techniques  for  obtaining  representative  samples  of  debris  from  hazardous-waste  sites.  As  stated  and 
discussed  by  Rupp  (1990),  the  problems  include  1)  obtaining  a  representative  sample  from  a  mix  of 
materials  of  various  sizes  and  compositions,  2)  characterizing  the  contamination  of  large  items  in  a  way 
that  has  meaning  for  a  health  risk  assessment,  and  3)  subsampling  from  mixtures  of  large  objects  to 
produce  the  small-volume  samples  required  by  the  analytical  laboratory.  The  basic  question  is  whether 
a  defined  unit  of  material,  e.g.,  a  barrel,  contains  areas  of  contamination  that  exceed  action  levels.  A 
related  problem  is  how  to  reduce  the  number  of  barrels  that  need  to  be  opened  and  characterized.  A 
sampling  approach,  such  as  acceptance  sampling  (Schilling,  1982),  might  be  feasible  to  resolve  the  latter 
problem.  But  generally,  the  solution  of  these  problems  should  be  a  team  effort,  wherein  DQOs  are 
established  first,  followed  by  studies  to  determine  sampling  and  inspection  approaches  that  meet  the 
DQOs. 

Coping  with  "Less-Than"  Data 

Frequently  the  concentration  of  a  contaminant  in  a  field  sample  cannot  be  quantified,  in  which 
case  it  may  be  reported  as  a  nondetectable  or  "less-than"  value.  When  such  less-than  values  are  present, 
it  is  difficult  to  obtain  valid  estimates  of  important  parameters,  such  as  mean  concentrations,  or  to  conduct 
valid  statistical  tests  to  identify  changes  in  concentrations  over  time  or  determine  compliance  with  clean-up 
standards.  The  common  practice  of  replacing  less-than  values  with  zeros  or  other  fabricated  values  can 
lead  to  highly  misleading  results.  To  avoid  this  problem,  it  is  necessary  to  use  special  statistical  methods 
(Helsel,  1990;  Gilbert,  1987).  The  need  for  additional  statistical  methods  for  such  cases  is  discussed  by 
Lambert,  et  al.  (1991)  who  also  introduce  new  tools  (tiie  "probabihty  of  acceptance"  and  the  "probability 
of  detection")  to  describe  which  measurements  and  field  concentrations  are  detectable. 
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Using  Uncertainty  and  Sensitivity  Analyses 

The  assessment  of  risks  associated  with  various  clean-up  scenarios  and  technologies  requires  the 
use  of  environmental  transport  and  risk  models,  the  predictions  of  which  are  often  highly  uncertain.  Yet 
this  uncertainty  may  not  be  explicitly  taken  into  account  when  formulating  DQOs  and  making  decisions. 
User-friendly  computer  codes  are  available  for  conducting  Monte  Carlo  uncertainty  analyses  (Iman  and 
Shortencarier,  1984)  and  sensitivity  analyses  (Iman,  et  al.  1985)  to  quantify  the  uncertainty  in  model 
predictions  and  to  identify  model  parameters  that  have  a  big  impact  on  model  predictions.  Use  of  these 
methods  may  be  considered  if  DQOs  require  precise  estimates  of  uncertainty  and  when  decisions  must  be 
made  about  which  model  parameters  should  be  refined  to  reduce  uncertainty.  The  use  of  uncertainty  and 
sensitivity  analyses  has  the  additional  benefit  of  clarifying  the  sources  of  uncertainty  in  the  model  and 
model  parameters.  This  identification  process  also  helps  develop  a  more  thorough  understanding  of  the 
uncertainties  present  in  the  system  and  where  they  reside.  IAEA  (1989)  has  described  procedures  for 
evaluating  the  reliability  of  predictions  made  by  environmental  transport  models,  including  uncertainty  and 
sensitivity  analyses  and  validation  studies.  Finkel  (1990)  has  discussed  uncertainty  in  risk  management 
for  decision-makers. 
Using  Nonparametric  Statistical  Tests 

The  standard  assumption  that  underlies  many  statistical  tests  of  hypotheses  is  that  the  data  are 
normally  (Gaussian)  distributed,  which  is,  however,  not  usually  the  case  with  waste-site  data.  In  such 
situations,  nonparametric  statistical  tests  should  be  considered.  For  example,  consider  the  problem  of 
testing  for  attainment  of  risk-based  or  background-based  soil  clean-up  standards  at  a  remediated  waste  site. 
The  nonparametric  Wilcoxon  Rank  Sum  test  and  Quantile  test  (Gilbert  and  Simpson,  1990)  can  be  used 
for  the  background-based  case.  For  the  risk-based  case,  nonparametric  tests  for  proportions  based  on  the 
binomial  distribution  can  be  used  (EPA,  1989c).  Nonparametric  tests  are  not  a  cure-all.  For  example,  the 
data  must  be  uncorrelated  in  order  for  the  nonparametric  test  results  to  be  valid,  the  same  requirement  as 
for  standard  parametric  tests.  However,  nonparametric  tests  can  be  more  powerful  (i.e.,  have  a  smaller 
false-negative  error  rate)  than  parametric  tests  when  the  assumption  of  normality  is  not  valid.  Also,  some 
nonparametric  tests  can  be  used  even  when  a  moderate  number  of  less-than  values  are  present  (Gilbert, 
1987),  if  all  less-than  values  are  less  than  the  smallest  detected  value.  The  Quantile  test  can  be  used  even 
when  a  large  proportion  of  the  data  are  less-than  values.  Regulators  and  the  planning  teams  for  ER  need 
to  become  familiar  with  nonparametric  tests  and  their  advantages  and  disadvantages. 
Communicating  with  the  Public 

The  public  supports  ER  with  tax  dollars.  Every  effort  should  therefore  be  made  to  effectively 
communicate  the  plans,  methods,  results,  and  implications  of  results  to  the  public  in  forms  that  the  average 
person  can  understand.  This  is  a  formidable  challenge.  Statistical  graphics  and  geographical  information 
systems  (GIS)  are  tools  that  can  contribute  to  this  communication  effort  (Tzemos,  et  al.  1991;  Dangermond 
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and  Hamden,  1990).  For  example,  a  GIS  can  be  used  to  integrate  hazardous- waste  site  data  with 
associated  geographic  information  via  map  overlays.  There  is  also  potential  for  integrating  statistical 
analyses  and  modehng  of  spatial  data  into  GIS  software,  although  present  capabilities  are  limited  (Bailey, 
1990).  This  is  an  area  where  more  work  is  needed  to  assess  what  is  required  and  what  the  costs  might 
be  to  develop  the  methodology.  Also,  statisticians  can  take  mysticism  out  of  statistical  techniques  by 
avoiding  jargon,  going  back  to  first  principles,  and  using  intuitive  descriptions  and  examples. 

Conclusions 

The  problems  of  cost  and  uncertainty  associated  with  ER  can  be  tackled  using  data  quality 
objectives,  the  observational  approach,  and  statistical  designs  and  analyses  developed  by  the  statistician 
in  collaboration  with  other  members  of  the  ER  planning  team.  A  number  of  statistical  tools  are  currently 
available,  but  the  development  of  additional  tools  is  needed.  This  paper  has  described  a  few  of  the  issues 
and  needs  related  to  the  application  of  statistics  to  site  characterization  and  ER  at  hazardous  waste  sites. 
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Some  Problems  in  Statistical  Consulting 
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Abstract 

The  role  of  consulting  statisticians  in  agricultural  research  has  changed  in  recent  years,  with 
collaboration  and  advising  on  research  studies  now  being  their  major  activity.  The  trend  for  researchers 
to  undertake  some  or  all  of  their  own  statistical  work  has  merit  if  they  have  a  good  knowledge  of 
statistical  methods  and  analysis  techniques,  but  can  lead  to  lower  research  standards  otherwise.  Some 
guidelines  for  the  handling  of  statistical  problems  are  discussed. 

Introduction 

Corresponding  to  advances  that  have  taken  place  in  computer  technology  in  recent  years,  there 
has  been  an  increase  in  the  number  of  users  of  statistical  methods  in  agricultural  research.  These  include 
people  with  considerable  formal  training  and  experience  in  statistics,  as  well  as  those  whose  main  expertise 
is  in  another  field,  but  within  which  statistics  plays  an  important  part.  Before  the  arrival  of 
microcomputers,  the  bulk  of  statistical  work  in  agricultural  research  was  handled  using  large  mainframe 
computer  systems,  and  often  through  a  professional  statistician  or  individual  who  had  sufficient  statistical 
knowledge  to  be  designated  as  the  "local  statistician". 

Today,  with  the  proliferation  of  computer  terminals,  microcomputers  and  statistical  software, 
increased  amounts  of  information  to  be  analyzed,  and  limited  numbers  and  availability  of  consulting 
stafisticians,  more  people  with  varying  degrees  of  statistical  training  are  carrying  out  their  own  statistical 
work.  This  has  led  to  a  change  in  the  role  of  the  statistician,  and  the  real  danger  of  misapplication  of 
statistical  methods  by  researchers  with  limited  statistical  backgrounds.  I  would  like  to  provide  an 
overview  of  statistical  consulting  at  an  agricultural  research  station  and  indicate  some  general  guidelines 
for  and  problems  encountered  when  handling  statistical  problems. 

Statistics  at  the  Lethbridge  Research  Station 

The  Lethbridge  Research  Station  is  an  agricultural  research  station  in  Agriculmre  Canada  that  has 
research  programs  in  Crop  Science,  Livestock,  and  Soil  Science  sections.  About  65  scientists  carry  out 
research  studies,  and  statistical  considerations  play  an  important  part  in  many  of  these.  The  Station  has 
3  statisticians  and  a  data  analyst  who  work  with  the  scientists  on  statistical  aspects  of  their  studies.  Most 
of  the  statistical  processing  is  carried  out  using  SAS  (SAS  Institute,  Inc.)  software  on  a  VAX  system,  but 
other  statistical  software,  on  the  VAX  and  microcomputers,  is  also  used. 
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The  statisticians  have  limited  to  major  input  in  statistical  aspects  of  specific  research  studies.  The 
statistician's  role  can  be  as  a  collaborator,  advisor,  or  to  provide  a  service.  As  a  collaborator,  the 
statistician  plays  an  essential  part  in  the  study,  working  closely  with  the  scientist  and  carrying  out  most 
of  the  statistical  work.  As  an  advisor,  tiie  statistician  is  consulted  about  various  aspects  of  the  design  and 
analysis,  but  the  major  burden  of  the  work  is  done  by  the  scientist.  When  providing  a  service,  entire 
routine  tasks  are  taken  on  by  the  statistician,  and  there  is  a  low  level  of  personal  involvement  in  the  study. 
In  recent  years,  the  service  role  has  diminished  to  a  low  level,  with  researchers  performing  routine 
analyses  and  the  consulting  statistician  spending  most  of  his  time  as  a  collaborator  or  advisor.  Genuine 
collaboration,  where  the  statistical  problem  is  challenging  and  the  scientist  and  statistician  pool  their  talents 
and  expertise  to  produce  high  quality  research,  is  the  most  rewarding  role  for  the  consulting  statistician. 

Although  it  is  beneficial  for  researchers  to  use  statistical  methods,  this  can  result  in  inadequate 
analyses,  misapplication  of  statistical  meUiods,  and  inefficient  use  of  computer  resources  if  they  do  not 
have  a  good  knowledge  of  statistical  metiiods  and  analysis  techniques.  Researchers  need  to  use  statistical 
metiiods  frequentiy  to  retain  tiieir  sidlls  and  develop  the  ability  to  recognize  statistical  problems  and 
processes  that  they  can  handle  themselves,  and  those  for  which  professional  assistance  is  required.  In 
addition,  the  researcher  must  consider  the  time  requirement  for  statistical  analysis  so  that  it  does  not  take 
too  much  time  away  from  his  own  discipline.  Easy  to  read  documentation  on  the  use  of  statistical 
software  combined  with  one-to-one  assistance  from  statistical  personnel  is  helpful  for  improving  their 
effectiveness  and  self-interest.  Even  when  advice  is  obtained  from  a  professional  statistician,  inappropriate 
analyses  can  occur  if  the  statistician  does  not  have  sufficient  information  on  or  time  to  digest  a  problem, 
often  as  a  result  of  work  pressures  and  a  backlog  of  requests. 

Some  guidelines  for  handling  statistical  problems 

When  working  with  scientists  on  a  research  study,  the  statistician  can  enter  the  study  at  the 
planning,  analysis,  or  interpretation,  presentation,  and  publication  stages.  Although  it  is  ideal  to  enter 
during  the  planning  and  experimental  design  stage,  and  be  fully  involved  at  all  subsequent  stages,  the 
statistician  often  enters  at  the  analysis  or  later  stages.  It  is  extremely  important  that  the  statistician  discuss 
tiie  problem  in  detail  with  the  scientist  to  get  an  accurate  and  complete  description  and  understanding  of 
those  parts  of  the  research  that  have  an  implication  on  the  design  and  statistical  analysis.  The  objectives 
should  be  clearly  formulated,  and  plans  for  the  design  and  analysis  outiined.  Good  dociunentation  of  key 
discussions  with  the  scientist  at  all  stages  of  the  study  are  valuable  for  the  efficient  handling  of  problems. 
For  complex  projects,  a  written  summary  of  the  statistician's  understanding  of  the  project  and  the  agreed 
on  course  of  action  should  be  sent  to  the  scientist. 

Prior  to  doing  a  detailed  statistical  analysis  of  the  data,  considerable  time  should  be  spent  careftilly 
scrutinizing  the  data,  both  visually  and  using  computer  summaries  and  plots,  to  check  for  entry  errors, 
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discrepant  values,  heterogenous  variance,  non-linearity,  and  so  on.  Initial  analyses  on  the  data  and  another 
review  of  the  objectives  with  the  scientist  will  clarify  the  analyses  that  are  really  important  and  the 
approach  to  take  for  obtaining  these.  Detailed  consideration  has  to  be  given  as  to  whether  data 
transformations  are  necessary,  the  choice  of  statistical  techniques  to  be  used  when  there  are  different 
possibilities,  and  whether  standard  or  specialized  statistical  methods  are  needed. 

Commonly  used  statistical  methods  found  in  standard  texts  like  Bliss  (1967),  Cochran  and  Cox 
(1957),  Draper  and  Smith  (1981),  Heiss  (1981),  Milliken  and  Johnson  (1984),  Snedecor  and  Cochran 
(1980),  Steel  and  Torrie  (1980)  and  Winer  (1962)  can  be  used  in  the  majority  of  research  problems  and 
are  more  likely  to  be  understood  by  researchers.  However,  there  are  times  when  nonstandard  or  highly 
complex  statistical  procedures  are  needed  to  meet  the  objectives.  For  example,  at  an  agricultural  research 
station  plant  breeding  studies  may  be  carried  out  to  evaluate  many  varieties  and  lines  for  various  characters 
such  as  yield  and  plant  height.  It  is  often  of  interest  to  determine  groupings  of  the  varieties  and  lines  that 
respond  similarly  at  different  locations,  over  years,  exhibit  a  similar  growth  response,  and  so  on.  A  joint 
regression  technique  (Finlay  and  Wilkinson  1963),  which  was  expanded  to  include  clustering  (Lin  and 
Thompson  1975),  is  useful  for  investigating  complicated  two  factor  interaction  structures  such  as  these. 

Once  the  statistical  analyses  have  been  obtained,  the  output  should  always  be  checked  very 
carefully  to  ensure  that  the  analyses  are  correct  and  that  the  statistics  produced  correspond  to  what  the  raw 
data  indicate.  Statistical  software  packages  can  give  incorrect  or  inappropriate  results  that  only  become 
apparent  with  close  examination  of  the  output.  The  final  analyses  should  be  summarized  using  tables  and 
diagrams  that  are  fully  informative  of  the  data  from  which  they  have  been  derived,  and  careftilly 
interpreted  with  the  original  objectives  in  mind.  A  written  report  of  the  statistical  methods  used  and 
results  obtained  should  be  prepared  for  the  scientist  so  that  they  are  clearly  understood  and  properly 
reported.  Drafts  and  galley  proofs  of  the  publication  must  be  reviewed  carefully  to  ensure  that  what  is 
published  is  what  was  intended. 

Summary 

The  ways  in  which  statistical  design  and  analysis  problems  are  handled  has  changed  in  recent 
years.  With  the  increased  availability  and  use  of  computers,  consulting  statisticians  are  assuming  primarily 
collaborative  or  advisory  roles,  and  more  researchers  are  carrying  out  their  own  statistical  work.  Although 
this  has  merit,  the  use  of  statistical  methods  by  individuals  with  a  limited  knowledge  of  statistical  methods 
and  analysis  techniques  can  lead  to  insufficient  or  incorrect  analyses,  and  important  and  interesting  features 
of  data  may  not  be  brought  out.  It  is  important  for  researchers  to  develop  skills  to  identify  problems 
where  the  appropriate  statistical  analysis  is  outside  their  knowledge  or  ability.  Some  guidelines  for  the 
handling  statistical  problems  are  discussed. 
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SAS  for  OS/2 


Serge  Dupuis 
Software  Support 
Alberta  Public  Works  Supplies  and  Services 
6950-113  Street 
Edmonton,  Alberta  T6H  5T6 

Abstract 

The  Statistical  Analysis  system  (SAS)  software  has  been  developed  to  operate  on  several  different 
computers.  This  includes  Mainframes,  Workstations  and  Personal  Computers,  each  with  its  own  operating 
system  (TSO,  VAX,  DOS,  UNIX  and  OS/2).  This  provides  a  great  deal  of  flexibility.  Users  can  run 
small  analyses  on  DOS,  access  large  databases  on  mainframes  and  move  their  SAS  Programs  to  UNIX 
or  OS/2  for  large  statistical  or  graphics  applications.  Each  version  of  SAS  has  the  same  user  language 
and  procedures  and  is  different  only  in  the  exploitation  of  each  different  operating  system.  This  paper 
introduces  some  of  these  features. 

OS/2  features 

SAS  for  OS/2  is  extremely  flexible  by  exploiting  features  such  as  Named  pipes.  Dynamic  data 
exchange  (DDE)  and  others.  Many  of  these  are  detailed  in  a  recent  paper  in  the  IBM  Personal  Systems 
Developer  (Goldstein,  1991) 

Named  Pipes  allow  SAS  to  establish  many  types  of  communication  links.  Within  ONE  machine, 
multiple  SAS  sessions  can  be  running,  with  actions  of  one  SAS  session  controlling  another.  With  Named 
pipes,  custom  C  applications  can  link  to  SAS.  The  support  for  Named  Pipes  allows  for  truly  distributed 
processing  in  a  networked  environment.  For  example,  SAS  for  OS/2  on  one  machine  could  be  collecting 
real-time  data  from  a  serial  data  collection  device,  while  another  computer  on  the  network  could  receive 
and  analyse  the  data.  Named  pipes  support  is  built  into  the  SAS  system  and  is  independent  of  LAN 
software.  It  will  run  with  any  LAN  software  supporting  named  pipes,  such  as  NOVELL  for  OS/2,  IBM 
LAN  server  and  Microsoft  LAN  manager. 

Dynamic  Data  Exchange  (DDE)  is  useful  for  establishing  communications  with  other  applications 
supporting  DDE.  Hot  links  can  be  established  between  EXCELL  for  example  and  SAS.  As  data  is 
changed  in  the  spreadsheet,  SAS  can  analyse  it  directly. 

Memory  management 

SAS  is  designed  to  take  advantage  of  all  available  extended  memory,  up  to  1  gigabyte  of  virtual 
memory.  This  memory  is  used  for  both  executable  code  and  data.  The  code  generator  used  in  the  DATA 
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step  and  the  PROC  code  have  been  optimized  for  OS/2.  Large  statistical  and  graphics  applications  will 
show  a  marked  improvement  in  performance  with  this  increased  memory. 

Multitasking 

OS/2  offers  full  multitasking,  not  time  sliced  tasking  as  used  in  DOS/WINDOWS.  This  means  that 
any  one  application  does  not  'take  over'  the  processor,  as  it  does  under  WINDOWS.  Also,  there  is  full 
protection  between  applications,  so  that  if  one  application  fails,  the  entire  computer  does  not  need  to  be 
'booted'. 

Disk  Access 

OS/2  can  use  the  High  Performance  File  System  (HPFS)  which  is  optimized  for  multitasking  and 
input/output  throughput.  SAS  under  OS/2  shows  25  to  50%  increase  in  performance  using  HPFS,  over 
the  standard  DOS  File  Allocation  Table  standard. 

Comparison  with  Windows  3.0 

S AS/PC  under  Microsoft  Windows  3.0  will  be  available  over  the  next  year,  however  OS/2  and 
Unix  are  likely  to  be  preferred  by  high-end  users  for  complex  applications.  A  paper  comparing  OS/2  and 
Windows  3.0  details  some  important  differences  in  memory  management,  multitasking  and 
communications  (Gates,  1990).  Under  Windows,  data  and  allocated  memory  are  shared  amongst 
applications  which  can  lead  to  data  corruption.  Under  OS/2,  there  can  be  up  to  two  or  three  SAS  sessions, 
while  under  Windows,  only  one  SAS  session  is  recommended.  Limits  to  SAS  datasets  are  the  same  for 
both  Windows  and  OS/2,  there  is  no  limit  to  the  number  of  observations  but  each  observation  has  a 
maximum  of  30,700  variables  of  2  bytes  each  (or  less  variables  of  greater  length). 
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Analysis  of  Air  Monitoring  Data 


L.Z.  Florence  and  H.  Bertram 
Alberta  Environmental  Centre 

Postal  Bag  4000 
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Abstract 

Pollution  monitoring  has  resulted  in  many  measurement  programs.  Monitoring  networks  are  used 
mmnly  to  understand  dispersion  processes  and  measure  deposition  rates  and/or  levels.  This  paper  is 
prepared  as  an  introduction  to  the  types  of  networks,  specifically  in  Alberta,  and  the  characteristics  of 
spatial  and  temporal  data  which  one  might  expect  to  analyze  or  sunmiarize;  examples  come  from  dry 
sulphate  and  SO2  deposition  networks.  An  outline  of  methods  and  literature  sources  is  presented  for  the 
statistical  treatment  of  air  monitoring  data. 

Introduction 

Public  concern  for  air  quality  and  atmospheric  deposition  has  resulted  in  numerous  measurement 
programs,  ranging  from  single  station  short-term  research  projects  to  long-term  nation-wide  monitoring 
programs.  Current  technology  limits  the  scope  of  the  monitoring  programs  to  determining  atmospheric 
concentrations;  however,  deposition  models  and  deposition  measurement  procedures  are  rapidly 
developing. 

The  monitoring  networks  are  used  for  at  least  two  purposes:  first,  to  help  in  our  understanding 
of  the  process  (emission,  transport  and  deposition),  and  second,  to  facilitate  a  specific  job  such  as  to  verify 
a  predicted  concentration  gradient,  or  deposition  rate.  Short  range  monitoring  is  usually  more  successful 
because  the  pollutant  does  not  have  time  to  undergo  chemical  transformation  or  deposition  to  the  earth's 
surface.  Long  range  monitoring  has  to  take  into  account  surface  deposition  and  both  primary  emissions 
and  secondary  emission  products,  such  as  the  formation  of  secondary  sulphate  in  the  form  of  ammonium 
sulphate  aerosol  from  sulphur  dioxide  emissions. 

This  paper  presents  discussion  and  examples  of  (1)  types  of  monitoring  networks,  (2) 
characteristics  of  monitoring  data,  and  (3)  methods  for  analyzing  monitoring  data. 

Types  of  Monitoring  Networks 

Monitoring  networks  tend  to  fall  into  three  broad  categories  (Whelpdale  1983):  exploratory,  where 
they  are  alerting  the  public  and  scientific  communities  to  the  problem  (ozone/nox/hydrocarbon  and 
greenhouse  gas  networks),  assessment,  where  methods  are  being  refined  and  the  extent  of  the  problem  is 
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being  defined  (such  as  the  acid  deposition  network)  and  regulatory,  where  the  effectiveness  of  the 
regulations  are  being  assessed  (such  as  the  urban  lead  network). 

Air  monitoring  data  for  Alberta  are  collected  from  sites  around  the  province  with  the  intent  of 
determining  the  effects  of  population  growth  and  industrialization  on  the  environment.  Pollutants  of 
concern  are  detemiined  by  their  toxic/hazardous  properties  and  the  ratio  of  anthropogenic  emissions  to 
natural  sources  of  emissions.  Industries,  and  the  provincial  government  as  a  quality  control  check, 
maintain  monitoring  networks  in  the  vicinity  of  major  processing  plants  as  part  of  the  regulatory  and 
licensing  requirements. 

Network  data  are  characterized  by  spatial  distribution  (location),  time  (date  sampled),  and  intensity 
(atmospheric  concentration).  The  type  of  sampling  procedure  must  be  taken  into  account  in  judging  the 
appropriateness  and  reliability  of  the  data.  For  example,  continuous  analyzers  are  used  to  detennine 
violations  of  the  one  hour  exposure  limit,  and  monthly  integrated  samplers  are  used  for  assessing  long- 
term  emission  trends  and  wind-rose  type  of  spatial  distributions. 

Alberta  Envirormient  maintains  networks  in  urban  areas;  for  example,  Edmonton  has  a  twenty- 
seven  site  sulphation  network  and  Fort  McMurray  has  a  six  site  network  (Myrick  and  Asquin  1991).  The 
Department  also  has  sites  which  act  as  quality  control  checks  around  major  processing  plants. 

Monitoring  SO2  Levels  in  a  Network  Using  Passive,  Integrated  Samples 

An  established  type  of  sampling  device  used  in  monitoring  networks  is  the  sulphation  candle.  The 
sulphation  candle,  which  has  been  in  use  for  the  past  fifty  years,  consists  of  a  lead  dioxide  (peroxide) 
paste,  using  gum  tragacanth  as  adhesive,  painted  in  a  cylindrical  strip  around  the  middle  of  a  glass  jar. 
The  candle  is  exposed  in  the  field  for  one  month,  where  the  lead  dioxide  reacts  with  ambient  sulphur 
dioxide,  hydrogen  sulphide  and  mercaptans  to  form  lead  sulphate.  After  exposure  the  candles  are  returned 
to  the  laboratory  for  analysis. 

A  refinement  of  this  device  is  the  Huey  plate,  which  is  more  convenient,  less  expensive  and  has 
some  improvement  in  accuracy.  Accuracy  is  improved  by  using  a  precise  amount  of  lead  dioxide  for  each 
sampler,  thus  eliminating  some  of  the  "blank"  variability.  Huey  plates  are  prepared  by  pipetting  a  paste 
consisting  of  lead  dioxide,  gum  tragacanth  and  cellulose  fibres  into  a  petri  dish  and  then  drying  the  plate 
in  an  oven.  The  sampler  is  then  sent  to  the  field,  and  exposed  open  side  down  in  a  convenient  holder. 

A  more  recent  development  in  passive  sulphur  dioxide  monitors  is  the  use  of  potassium 
carbonate/glycerol  on  cellulose  filter  paper  as  the  absorbent.  As  well,  diffusion  screens  are  now  being 
used  to  eliminate  wind  velocity  effects.  Extensive  networks  of  these  low  cost  devices  have  been 
established  in  Alberta.  The  data  gathered  from  these  devices  are  used  to  delineate  major  areas  of  sulphur 
dioxide  impact,  and  to  identify  trends  in  deposition  levels. 
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Analyzing  and  Summarizing  Monitoring  Data:  examples  using  sulphur  dioxide  (SO2)  emissions  data 
from  Alberta 

Gilbert  (1987),  and  references  therein,  discusses  statistical  methods  for  assessing  spatial  and 
temporal  distributions  of  monitoring  data.  The  reader  is  referred  to  these  sources  for  more  detailed 
discussion  of  the  methods  suggested  here. 
Spatial  Patterns 

Integrated  samplers  (total  monthly  deposition)  in  a  network  provide  data  that  can  be  summarized 
monthly  or  yearly.  Commonly  applied  methods  for  assessing  spatial  patterns  are  trend  surface  analysis 
and  kriging,  a  weighting  technique  utilizing  the  correlation  structure  among  locations  in  a  grid  or  network. 
An  example  of  how  average  annual  dry  sulphate  deposition  may  be  summarized  using  data  from  local 
networks  represented  by  concentration  isopleths  is  shown  in  Figure  1  (Fort  McMurrav  Regional  Air 
Quality  Task  Force,  1989)  for  the  Suncor  (N=40  sampling  sites)  and  Syncrude  (n=40  sampling  sites). 
Note  the  concentration  gradients  with  distance  from  each  plant  installation. 
Spatial  and  Temporal  Patterns 

Figure  2  shows  comparisons  among  four  of  the  27  SO2  sampling  locations  in  the  Edmonton 
network.  Each  observation  is  an  integrated,  monthly  concentration  of  SO2  from  January  1980  through 
December  1989.  These  plots,  therefore,  display  both  historical  temporal  variation  and  among  location 
(within  network)  variability. 

To  model  these  relatively  long  data  series,  a  Box- Jenkins  (Box  and  Jenkins  1976)  approach  could 
be  applied  (requiring  at  least  50  data  points),  using  software  such  as  SAS/ETS  (SAS  Instimte,  1988),  for 
detecting  trends  or  forecasting  future  deposition  rates. 

Referring  to  Figure  2,  we  might  conclude  that  all  four  locations  experienced  a  declining  trend  in 
SO2,  starting  with  the  1980  data  continuing  until  at  least  tiie  mid-1980's.  The  latter  half  of  the  1980's  data 
appear  to  fluctuate  in  a  random,  stationary  manner.  Gilbert  (1987),  chapters  16  and  17,  provide  a  very 
useful  guide  for  methodically  assessing  trends  in  monitoring  data.  Most  of  the  techniques  discussed  by 
him  are  robust  enough  that  outliers,  missing  and/or  non-detectable  data,  are  permissable. 

A  protocol  for  analyzing  the  data  shown  in  Figure  2  might  follow  these  steps: 

1.  First,  plot  tiie  data;  graphical  analysis  can  be  very  powerful  for  not  only  assessing  whether  data 

contain  patterns,  but  also  as  a  quality  control  for  detecting  data  errors  or  exceptional  observations 

which  should  receive  special  attention. 
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*    SUNCO«  SULPHATlON  STATION 
»    5YNCRU06  SUIPMATON  STATION 

Figure  1.  Estimates  of  average  annual  dry  sulphate  deposition  (kg  ha'^)  derived  from  information 
obtained  from  Suncor  and  Syncrude  sulphates  data  collected  from  1979  to  1983. 
(Figure  2,  p.  10,  Dabbs,  1985). 
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Figure  2.   S02  Concentration  Over  Time  At  Four  Edmonton  Locations 
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Figure  2.         SO2  Concentration  Over  Time  at  Four  Edmonton  Locations. 
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2.  Once  plotted,  for  example,  in  Figure  2,  we  may  wish  to  test  the  null  hypothesis  that  no  trends 
exist  over  time  and  among  locations.  Total,  global  variation  in  SO2  can  be  partitioned  and  the  null 
hypothesis  tested,  as  seen  in  Table  1,  Part  A.  Because  we  would  expect  only  a  0.021  chance  of 
obtaining  a  larger  chi-square  value  for  "Station",  we  can  ignore  "Trend"  and  conclude  that 
statistically  significant  different  trend  directions  occurred  between  at  least  two  of  the  four 
Edmonton  stations  during  the  10  year  period  plotted  in  Figure  2.  Table  1,  Part  B,  indicates  that 
all  four  locations  experienced  significant  trends  and,  as  might  be  expected  from  Figure  2,  the 
weakest  statistically  was  at  #19. 


Table  1.         Tests  for  homogeneity  of  trends  among  four  Edmonton  locations;  data  are  displayed  in 
Figure  2. 


Part  A.  Tests  for  homogeneity  of  trends  among  four  Edmonton  locations. 

Source 

Chi-square 

df 

Total 

143.62 

48 

0.000 

Homogeneity 

55.58 

47 

0.183 

Season 

10.38 

11 

0.496 

Station 

9.75 

3 

0.021 

Station-Season 

35.45 

33 

0.354 

Trend 

88.04 

1 

0.000 

Part  B.  Tests  for  trend  at  each  Edmonton  Location. 

Location 

Chi-square 

df 

P>f 

4 

25.96 

1 

0.000 

6 

40.94 

1 

0.000 

11 

26.28 

1 

0.000 

19 

4.60 

1 

0.032 

^  Probability  of  a  larger  chi-square  value. 

3.  A  further  point  of  inquiry  would  likely  be  whether  the  trends  remained  significant  upon 
adjusting  for  seasonal  influence  and,  if  so,  what  direction  were  the  trends?  Table  2,  Part  A,  shows 
that  seasonally  adjusted  trends  existed  and  all  were  significant  at  the  0.03  level  or  less.  Table  2, 
Part  B,  contains  the  Seasonal-Kendall  slope  coefficients  and  95%  confidence  intervals  for  each 
location;  as  might  be  predicted  from  examining  Figure  2,  all  four  slopes  are  negative  and  the 
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weakest  trend  is  at  Station  19.  We  can  conclude  that  the  relative,  median  rate  of  change  in  SO2 
concentration  at  location  #19  was  approximately  one  third  that  observed  at  the  other  three  stations. 


Table  2.  Tests  for  seasonally  adjusted  SO2  trends  at  four  Edmonton  locations. 


Part  A. 

Tests  for  seasonally  adjusted  trends. 

Location                             Seasonal  Kendall  Score 

P>  SK 

1 

4 

-5.08 

0.000 

6 

-6.38 

0.000 

11 

-5.10 

0.000 

19 

-2.13 

0.033 

Part  B. 

Slopes  and  their  confidence  limits. 

Location                                     Lower  Limit^ 

SK  Slope 

Upper  Limit 

4 

-0.011 

-0.008 

-0.005 

6 

-0.010 

-0.008 

-0.005 

11 

-0.013 

-0.009 

-0.006 

19 

-0.007 

-0.003 

0.000 

^  Probability  of  a  greater  Seasonal-Kendall  (SK)  score. 

^  Lower  and  Upper  95%  confidence  limits. 

An  obvious  question  arises  as  to  the  cause  for  the  decline  in  SO2  concentrations.  Two 
obvious  conclusions  are:  (i)  industry  was  successfully  reducing  emissions,  or  (ii)  due  to  the  severe 
downward  economic  change  in  Edmonton  about  1981-82,  emissions  declined  due  to  industry  shut- 
downs. Likely,  both  these  reasons  are  part  of  the  explanation,  but  a  third  is  also  very  important 
because  we  also  know  that  prior  to  1982,  a  different  lab  was  analyzing  the  SO2  samples,  and  that 
it  was  using  a  different  method  of  assay  than  was  used  after  the  change  in  labs.  This  type  of 
information  is  invaluable  when  attempting  to  attribute  cause-effect  in  trend  analysis. 

4.  Seasonal  trends  are  common  in  SO2  monitoring  data.  For  illustration,  10  means  and  medians, 
of  1980-89  monthly  data,  for  Edmonton  Location  #4,  are  shown  in  Figure  3  (smoothing  was  done 
with  a  cubic  spline  algorithm  available  in  SAS/STAT  [SAS  Institute,  1988]).  We  show  both  these 
measures  of  centrality  to  demonstrate  that  the  median  is  less  sensitive  to  extreme  values  than  the 
mean  (note  especially  August  through  November),  thus  better  using  all  data  rather  than  possibly 
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deleting  obvious  outliers.  When  the  data  are  symmetrically  distributed  (note  February,  April  and 
June)  the  mean  and  median  are  near  the  expected  values  for  normally  distributed  populations. 
Figure  3  clearly  shows  that  SO2  concentrations  increase  in  the  fall-winter  and  are  considerably 
reduced  in  the  spring-summer.  These  seasonal  trends  are  largely  due  to  factors  such  as  lower 
mixing  heights,  slower  rate  of  oxidation  to  sulphate  as  hydroxyl  radical  concentrations  decrease, 
and  lower  "washout"  during  the  winter  months,  than  during  summer,  resulting  in  greater  ambient 
SO2  deposition.  Similar  seasonality  was  observed  in  the  other  Edmonton  stations  (not  shown). 
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Figure  3.         Variability  in  SO2  concentration  about  the  mean  and  median,  over  10-year  period,  at  one 
sampling  location  4  in  Edmonton,  Alberta. 


Summary 

This  paper  has  been  an  introduction  to: 

1.  The  types  and  purposes  of  air  monitoring  networks; 

2.  Examples,  using  SO2  data  from  Alberta,  showing  how  data  may  be  distributed  over  time  (years 
and  seasons)  and  among  locations; 

3.  And  a  stepwise  approach,  albeit  not  comprehensive,  to  analyzing  monitoring  data. 
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